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We  present  three  numerical  models  of  accretion  from  radiation 
driven  stellar  winds  onto  compact  objects  in  massive  X-ray  binary 
systems.  The  wind  is  given  a  velocity  profile  consistent  with  a 
radiatively  driven  wind,  and  a  ^negative  mass"  gravitational 
potential  is  derived  from  this  profile  to  represent  the  wind  driving 
force  in  the  hydrodynamic  equations.  An  X-ray  heating  model  is  used 
which  determines  the  X-ray  heating  time  from  the  Compton  heating  time 
and  the  known  steady  state  energies  for  optically  thin  gas 
illuminated  by  X-rays.  This  allows  X-ray  heating  to  be  included  in 
the  hydrodynamic  equations.  The  X-ray  luminosity  is  held 
proportional  to  the  accretion  rate,  assuming  that  the  gravitational 
potential  energy  released  is  equivalent  to  10%  of  the  infalling 
rest-mass  energy.  A  two-dimensional  Eulerian  computer  code  is  used 
to  solve  the  equations  of  motion.  Model  estimates  of  the  ionization 
structure,  accretion  rates  and  flow  characteristics,  and  the  effects 
of  thermal  instabilities  are  discussed.  The  impact  of  the  X-ray 
radiation  on  the  wind  driving  force  is  demonstrated.  Results 
indicate  a  possible  mechanism  for  slow  X-ray  flares,  such  as  observed 
in  4U1 700-37. 

I\ 
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NUMERICAL  STUDIES  OF  GRAVITATIONAL  ACCRETION  FROM 
X-RAY  HEATED  STELLAR  WINDS 

I.  INTRODUCTION 

One  category  of  well  studied  cosmic  X-ray  sources  are  those 
found  in  binary  star  systems.  The  most  luminous  of  these  objects  are 
the  massive  X-ray  binary  stars.  These  massive  X-ray  binaries 
typically  contain  an  0  or  B  class  supergiant  primary  with  a  mass  of 
20  Mg  or  greater.  The  secondary  is  a  compact  object,  such  as  a  white 
dwarf,  neutron  star,  or,  possibly,  a  black  hole,  with  a  mass  of  10  Mg 
or  less.  The  secondary  is  the  source  of  the  X-ray  emission  with 

on  i 

luminosities  ranging  from  10  to  10  erg  sec"  (see  the  recent 

source  compilation  by  Amnuel,  Guseinov,  and  Rakhamimov  1979,  and 

references  cited  therein).  Observations  of  0  and  B  supergiants  have 

shown  that  these  stars  often  have  strong  stellar  winds  with  mass  loss 
-4  -1 

rates  of  up  to  10  Mg  yr  and  terminal  wind  velocities  on  the  order 
of  10  km  sec  (see,  for  example,  Weymann  1963;  Morton  1967a, b, 
1976;  Hutchings  1976,  1980;  Conti  1978a, b).  Such  observations  have 
lent  credence  to  the  generally  accepted  model  that  the  X-ray  emission 
from  these  binary  systems  is  powered  by  gravitational  accretion  of 
material  from  the  stellar  wind  onto  the  secondary.  This  model  was 
first  proposed  by  Davidson  and  Ostriker  (1973)  (see  also  Lamb, 
Pethick  and  Pines  1973). 

Development  of  our  current  understanding  of  gravitational 
accretion  began  with  the  work  of  Hoyle  and  Lyttleton  (1939).  They 
examined  the  gravitational  capture  of  material  by  a  body  moving 


supersonically  through  an  intergalactic  medium.  This  study  led  to 
the  crude,  but  useful,  analytic  line  accretion  model  (see  also  Bondi 
and  Hoyle  1944)  in  which  the  accretion  rate  was  estimated  using  the 
velocity  and  density  of  the  flow  and  the  mass  of  the  gravitating 
body.  Bondi  (1952)  considered  spherically  symmetric  accretion  in 
some  detail  under  the  assumption  that  the  gas  was  polytropic.  Hunt 
(1971,  1979)  generated  numerical  models  which  examined  gravitational 
accretion  near  the  transition  region  between  subsonic  and  supersonic 
flows  of  an  adiabatic  gas.  Hoffman  (1979)  created  a  non-radiative 
hydrodynamical  model  of  accretion  flow  in  an  isothermal  gas. 
Isothermal  flow  is  much  closer  to  the  problem  of  an  X-ray  source 
accreting  material.  The  X-rays  produce  signifcant  Compton  heating 
and  cooling  effects  near  the  secondary.  Heating  of  the  gas  occurs 
faster  than  the  material  flow  times,  negating  the  assumptions  of  an 
adiabatic  gas.  Hoffman  found  a  dense,  fluctuating  wake  downstream 
from  the  body.  Fluctuations  also  appeared  in  the  accreting  column  of 
matter  falling  back  to  the  secondary.  This  resulted  in  corresponding 
fluctuations  in  the  expected  X-ray  luminosity. 

Numerous  authors  (Mestel  1954;  Shvartsman  1971;  Buff  and  McCray 
1974a)  have  examined  cases  of  ideal  spherical  accretion  under  the 
influence  of  radiation.  They  found  that  the  X-ray  radiation  could 
seriously  affect  the  accretion  flow  under  the  combined  effects  of 
radiation  pressure  and  the  increases  in  gas  pressure  due  to  heating 
of  the  gas.  These  increased  pressures  led  to  severe  reductions  in 
the  amount  of  material  accreted.  Carlberg  (1978),  testing  radiation 


effects  analytically  in  the  line  accretion  model,  found  that  a 
turbulent  wake  was  probable,  and  that  strong  perturbations  of  the 
accretion  flow  were  possible  depending  on  the  physical  parameters  of 
the  flow.  He  concluded  that  hydrodynamic  calculations  would  be 
necessary  to  determine  the  actual  physical  states  of  the  accretion 
flow.  Cowie,  Ostriker,  and  Stark  (1978)  performed  a  time  dependent 
calculation  of  spherically  symmetric  accretion  onto  compact  X-ray 
objects.  They  demonstrated  strong  modulation  of  the  accretion  flow 
due  to  X-ray  heating  which  led  to  flaring  of  the  X-ray  luminosity. 
They  concluded  that  their  model  was  of  limited  application  to  X-ray 
binaries  because  of  the  strong  asymmetries  introduced  by  the  high 
wind  and  orbital  velocities. 

The  effects  of  X-ray  heating  and  ionization  of  diffuse  gases 
were  investigated  by  Tarter,  Tucker  and  Salpeter  (1969).  More 
recently.  Buff  and  McCray  (1974a)  and  Hatchett,  Buff,  and  McCray 
(1976)  investigated  the  heating  and  ionization  of  diffuse  gas  by 
X-rays  including  the  effects  of  Compton  heating.  They  pointed  out 
possible  multi-valuedness  of  the  temperature  function  and  the 
possible  formation  of  thermal  instabilities  (Field  1965).  In 
particular,  McCray  and  Hatchett  (1975)  considered  the  effects  of 
X-ray  heating  on  a  wind  and  deduced  that  such  flows  could  have 
multiple  pressure-density  states.  Alme  and  Wilson  (1975,  1976)  have 
demonstrated  that  instabilities  will  occur  in  X-ray  driven  flows. 
Hatchett  and  McCray  (1977)  examined  the  ionization  structure  in  the 
stellar  wind  of  an  X-ray  binary  assuming  an  unperturbed  wind  and 
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constant  luminosity. 

The  concept  of  stellar  winds  driven  by  radiation  pressure  was 
first  formulated  by  Holzer  and  Axford  (1970)  and  Lucy  and  Solomon 
(1970).  In  1975,  Castor,  Abbott  and  Klein  produced  a  detailed  model 
for  a  stellar  wind  which  included  contributions  to  the  driving  force 
by  line  absorption  in  trace  elements  (for  a  review,  see  Casinelli 
1979).  They  found  that  the  principle  contributions  to  the  force  were 
due  to  L-shell  absorption  lines  in  the  CNO  series  of  elements.  They 
also  found  a  velocity  law  for  the  wind  which  agreed  with  the  earlier 
studies.  But  most  significantly,  they  were  able  to  account  for  the 
wind  velocities  and  huge  mass  loss  rates  inferred  from  the 
observations  of  OB  supergiants.  The  success  of  this  concept  raises 
the  possibility  that  severe  interference  may  exist  between  the  wind 
and  the  X-ray  radiation  through  the  ionization  of  the  elements 
responsible  for  the  wind  driving  force. 

Since  the  accretion  model  of  powering  binary  X-ray  sources  is  so 
widely  accepted,  we  have  attempted  to  develop  a  model  which 
simultaneously  treats  the  hydrodyna  ,c  and  radiation  effects.  The 
model  is  intended  to  be  self-consistent  to  the  extent  that  the  X-ray 
luminosity  depends  directly  on  the  accretion  rate.  The  powering 
mechanism  assumes  gravitational  accretion  from  the  stellar  wind. 

In  this  study,  we  present  three  self-consistent  models.  By 
adopting  the  stellar  wind  model  presented  by  Castor,  Abbott,  and 
Klein  (1975)  we  are  able  to  include  the  force  driving  the  wind  as  a 
"negative  mass"  gravitational  potential  added  to  the  gravitational 
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potential  of  the  system.  This  wind  model  enables  us  to  account  for 
the  impact  of  the  X-ray  ionization  on  the  wind.  Hatchett,  Buff,  and 
McCray  (1976)  were  able  to  determine  the  ionization  equilibrium  for 
various  elements  in  a  gas  exposed  to  X-rays  as  a  function  of  a  single 
parameter.  We  make  use  of  the  parameter  in  our  model  to  determine 
the  ionization  state  of  the  elements  primarily  responsible  for  the 
wind  force.  The  ionization  state  in  turn  determines  the 
effectiveness  of  the  driving  force  included  in  the  hydrodynamic 
equations.  The  X-ray  luminosity  is  held  proportional  to  the 
accretion  rate,  thus  allowing  the  luminosity  to  vary  with  any 
disturbances  set  up  in  the  flow  by  the  X-rays.  X-ray  heating  rates 
are  explicitly  calculated  from  the  hydrodynamic  variables.  The 
resulting  calculations  provide  an  insight  into  the  flow 
characteristics  and  the  overall  ionization  structures  of  the  wind. 

In  Chapter  II,  we  summarize  the  observed  characteristics  of 
several  well  known  X-ray  binaries,  and  describe  the  two  broad  classes 
of  such  objects.  The  generally  accepted  powering  mechanism  converts 
the  gravitational  potential  energy  of  infalling  matter  into  thermal 
energy.  We  discuss  the  two  most  likely  means  of  providing  the 
necessary  matter,  namely  potential  overflow  and  gravitational 
accretion  from  a  stellar  wind,  with  emphasis  on  the  latter  method. 
The  limiting  effects  of  radiation  pressure  are  presented  in  a 
discussion  of  the  Eddington  luminosity.  We  then  present  a  physical 
argument  due  to  van  den  Heuval  (1975)  which  provides  a  basis  for  the 
existence  of  the  two  classes  of  X-ray  binaries. 


5 


The  methods  used  to  model  the  physical  processes  in  X-ray  binary 
systems  are  presented  in  Chapter  III.  Here  we  present  the  equations 
of  motion  with  two  terms  added  for  the  X-ray  heating  and  the  stellar 
wind  driving  force,  along  with  their  derivations.  The  X-ray  heating 
effects  are  handled  in  terms  of  an  X-ray  heating  time  computed  from 
the  Compton  heating  time  and  from  the  known  steady  state  energies  of 
a  gas  exposed  to  X-ray  radiation.  The  wind  force  term  is  derived 
from  the  stellar  wind  model  present  by  Castor,  Abbott,  and  Klein 
(1975).  We  show  that  their  wind  velocity  law  can  be  used  to  give  an 
effective  potential  for  the  force  acting  on  the  wind.  Finally,  we 
discuss  the  Lax-Wendroff  two-step  finite  differencing  method  which  is 
used  to  numerically  solve  the  equations  of  motion. 

In  Chapter  IV,  we  present  the  tools  required  to  develop  the 
physical  parameters  for  our  models.  Since  we  require  the  luminosity 
to  be  proportional  to  the  accretion  rate,  we  derive  relationships 
which  give  the  needed  accretion  rates  and  primary  mass  loss  rates  in 
terms  of  the  other  binary  parameters.  We  then  present  the  physical 
arguments  used  to  determine  the  acceptability  of  the  derived 
parameters. 

Three  different  models  were  considered  in  this  study,  using  two 
different  binary  systems  as  guides  for  parameter  selection.  These 
two  systems  were  4U0900-40  (Vela  X-l)  and  4U1700-37.  The  parameters 
for  these  systems  are  given  in  Table  2.1.  The  first  two  models  were 
designed  to  provide  reasonable  continuity  between  previous 
hydrodynamic  studies  involving  constant  material  flow  past  an 
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accreting  body  and  models  which  allow  both  variations  in  wind 
ionization  and  variable  material  velocities.  The  X-rays  are  allowed 
to  heat  the  wind  material,  but  no  attempt  is  made  to  alter  the  wind 
ionization  structure  and  thereby  modify  the  wind  driving  force. 

We  present  Model  1  in  Chapter  V.  This  model  is  based  on  the 
4U0900-40  system  and  model  parameters  are  very  close  to  those 
observed  in  the  actual  system.  We  find  that  the  accretion  rates 
obtained  are  a  factor  of  2  below  those  predicted  by  the  line 
accretion  model.  The  wind  driving  force  is  found  to  dominate  the 
secondary's  gravitational  force  beyond  a  certain  distance  from  the 
secondary.  We  also  demonstrate  that  the  flow  is  dominated  by 
hydrodynamics  rather  than  X-ray  heating  effects.  The  X-ray  heating 
is  shown  to  introduce  a  thermal  instability  into  the  flow.  It  also 
supresses  the  formation  of  an  accretion  column,  causing  most  of  the 
accreted  material  to  flow  in  from  the  side  regions  of  the  wake. 

Model  2  is  presented  in  Chapter  VI.  While  this  model  is  also 
based  on  the  4U0900-40  system,  the  model  parameters  were  adjusted  to 
place  it  in  a  state  such  that  the  flow  is  dominated  by  the  X-ray 
heating  rather  than  the  hydrodynamics.  The  resulting  model  shows  a 
large  deviation  from  the  behavior  suggested  by  the  line  accretion 
model.  The  accretion  rate  is  severely  impaired,  and  the  model 
developes  a  very  hot,  rarified  region  in  the  central  core  of  the 
wake.  The  model  also  has  a  thermal  instability  present  in  the  flow 
which  is  similiar  to  that  found  in  Model  1. 

Prior  to  generating  a  model  which  would  allow  the  wind  force  to 
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vary,  we  produce  a  simplified  ballistic  model  of  the  wind  flow  under 
the  impact  of  extreme  X-ray  ionization.  The  study,  presented  in 
Chapter  VII,  indicates  that  even  severe  modification  of  the  wind 
force  should  not  significantly  reduce  the  amount  of  material 
available  for  accretion  by  the  secondary.  It  does  demonstrate  that 
severe  three-dimensional  effects  are  present  in  the  flow  for  high 
X-ray  luminosities.  We  utilize  the  study  of  ionization  states  in  a 
wind  exposed  to  X-rays  presented  by  Hatchett,  Buff,  and  McCray  (1976) 
to  estimate  a  value  for  our  heating  parameter  at  which  the  wind  is 
considered  turned  off.  The  results  of  these  two  studies  provide 
additional  constraints  for  our  Model  3. 

In  Chapter  VIII  we  develop  the  parameters  for  Model  3  and 
present  the  results  of  our  calculation.  Due  to  the  restrictions  we 
develop  in  Chapter  VII,  we  select  4U1700-37  as  our  base  system.  The 
model  parameters  are  constrained  to  be  in  close  agreement  with  the 
observed  parameters  of  the  system.  Model  3  differs  from  the  previous 
two  in  that  we  allow  the  wind  force  to  be  affected  by  the  X-ray 
radiation.  The  ability  to  turn  off  the  wind  force  sets  up  a  feed 
back  mechanism  in  the  wind  flow.  This  mechanism  leads  to  large 
variations  in  the  wind  density.  As  these  density  fluctuations  pass 
by  the  secondary,  they  cause  peaks  in  the  accretion  rate  which  result 
in  slow  X-ray  flares.  The  flares  show  a  variation  of  3  orders  of 
magnitude  over  a  time  scale  of  about  2  hours.  This  flaring  can  be 
compared  with  the  slow  flares  observed  in  4U1700-37. 

Finally,  in  Chapter  IX,  we  discuss  the  results  of  our  model 
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calculations.  The  scaling  of  our  model  parameters  to  match  other 
systems  is  discussed.  We  perform  calculations  of  the  optical  depths 
in  the  models  and  find  that  the  models  become  optically  thick  in 
regions  far  from  the  secondary,  particularily  in  the  case  of  Model  3. 
We  then  discuss  the  weaknesses  of  Model  3,  and  compare  our  results  to 
some  of  the  most  recent  work  in  the  field.  In  particular,  we  compare 
our  models  to  those  produced  by  Livio,  et  al  (1979)  and  Hoffman 
(1979).  We  also  compare  our  results  to  Carlberg's  (1978)  analytic 
study  of  radiation  effects  in  the  line  accretion  model. 
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II.  THE  X-RAY  BINARY  STARS 
A.  Observational  Summary 

Following  Blumenthal  and  Tucker  (1974)  we  define  a  compact  X-ray 
source  as  one  identified  with  a  stellar  object,  or  one  which  exhibits 
variability  in  its  X-ray  luminosity  by  a  factor  of  2  or  more  on  a 
time  scale  of  days  or  less.  One  important  sub-group  is  composed  of 
the  X-ray  binary  stars.  In  these  systems,  one  of  the  stars  is  a 
normal  star  (hereafter  referred  to  as  the  primary)  and  the  other  (the 
secondary)  is  a  compact  object  such  as  a  white  dwarf,  a  neutron  star, 

OC  TO 

rr  a  black  hole.  They  typically  have  luminosities  of  10  to  10 
ergs  sec"*  in  the  2-10  keV  portion  of  the  X-ray  spectrum.  The 
characteristics  of  some  of  the  better  studied  of  these  stars  are 
$iven  in  Table  2.1.  Except  for  the  X-ray  pulsars  which  have  sharply 
defined  periodic  variations  in  their  luminosities,  we  see  that  the 
luminosities  of  these  objects  are  erratic  in  nature  showing  both  slow 
and  rapid  flaring  behavior  over  a  large  range  of  luminosities. 

By  examining  the  masses  of  the  objects  listed  in  Table  2.1,  we 
find  that  they  can  be  put  into  one  of  two  broad  classes.  One  class 
has  low  mass  primaries  of  around  2  Mg  or  less.  The  second  class 
contains  0  and  B  class  supergiant  primaries  with  masses  of  20  Mg  or 
more.  This  second  class  is  referred  to  as  the  massive  X-ray 
binaries.  In  both  classes,  the  secondary  masses  are  typically  1  or  2 
Mg.  With  the  recent  improvement  in  X-ray  observations,  there  is 
growing  evidence  for  a  third  class  of  X-ray  binary  (Amnuel,  Guseinov, 
and  Rakhamimov  1979).  This  class  has  luminosities  of  10^  erg  sec"* 
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approx.  0.5  keV  spectral  cut-off. 
Slow  flares  on  times  of  0.1  to  1  hr 
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or  less  and  have  indistinct,  extended,  infra-red  sources  as  their 
primaries. 

Since  the  compact  objects  are  the  sources  of  the  X-rays  in  these 
systems,  the  most  likely  method  of  powering  them  is  through  the 
release  of  gravitational  potential  energy  from  material  falling  onto 
the  object  (for  reviews  see  Blumenthal  and  Tucker  1974;  Giacconi 
1976;  Gurskey  1976).  The  amount  of  energy  released  in  the  form  of 
radiation  by  the  infalling  matter  was  shown  classically  by  Zel'dovich 
and  Shakura  (1969)  and  relati vistical ly  by  Alme  and  Wilson  (1973)  to 
be  on  the  order  of  10  to  20  percent  of  the  matter's  total 
rest-energy.  At  this  rate  of  energy  release,  we  need  only  provide  an 
accretion  rate  of  10“*^  to  10~®  Mg  yr~*  to  achieve  the  10^  to  10^ 
ergs  sec"*  luminosities  observed.  Such  mass  capture  rates  can  easily 
be  provided  by  one  of  two  different  mass  exchange  mechanisms  present 
in  binary  systems. 

B.  Mass  Exchange  Mechanisms 

In  1973,  Davidson  and  Qstriker  proposed  two  models  for  X-ray 
binaries  based  on  the  two  mass  exchange  mechanisms.  Tigure  2.1  shows 
schematically  the  processes  involved.  Figure  2.1(a)  depicts  the 
first  of  these  mechanisms,  that  of  potential  overflow.  The  figure-8 
shape  represents  the  critical  gravitational  potential  lobe,  which  is 
actually  a  surface  of  constant  gravitational  potential. 
Gravitational  potentials  at  greater  distances  than  the  critical 
potential  have  surfaces  which  surround  both  stars.  Potentials  closer 
than  the  critical  potential  have  surfaces  which  surround  each  star 
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Primary 


Figure  2.1(a).  Roche  Lobe  Overflow.  The  primary's  atmosphere 
(cross-hatching)  has  filled  the  critical  potential  lobe  and  is 
overflowing  onto  the  secondary.  (Not  shown  to  scale.) 


Figure  2.1(b).  Accretion  from  a  Stellar  Wind.  The  primary  does  not 
fill  the  critical  potential  lobe.  The  secondary  sweeps  up  stellar 
wind  material  from  its  orbit.  (Not  shown  to  scale.) 


individually.  To  express  these  surfaces  analytically,  let  us  take 
coordinates  x,  y,  and  z  such  that  the  z-axis  passes  through  the 
center  of  mass  of  the  system,  normal  to  the  orbital  plane,  and  the 
x-axis  lies  on  the  line-of-centers  with  the  neutron  star  on  the 
negative  x-axis.  The  gravitational  potential  at  any  point  is  then 
given  approximately  by  (Davidson  and  Ostriker  1973) 

$  =  -G  (Mx/rx  +  Mq^o  -  X(x,y)).  (2.1) 
Mx  and  r^  are  the  mass  of  the  secondary  and  the  distance  from  it. 
Likewise,  M  and  rQ  are  the  mass  of  the  primary  and  the  distance  from 
it.  X(x,y)  is  a  centrifugal  term  which  accounts  for  the  rotational 
characteristics  of  the  binary  system.  Let  us  consider  the  form  of 
X(x,y)  for  two  limiting  situations.  If  the  rotational  period  of  the 
primary  is  synchronized  to  the  secondary's  orbital  period,  we  are  in 
the  Roche  limit  and  X(x,y)  is  given  by 

X(x,y)  =  %  (Mq  +  Mx)(x2  +  y2)  a'3  ,  (2.2) 
where  a  is  the  binary  seperation.  The  second  limit  is  if  the  primary 
is  not  rotating  at  all.  This  case  is  called  the  tidal  limit,  and 
X(x,y)  is  given  by 

X(x,y)  =  x  Mx/a2.  (2.3) 
An  actual  binary  system  most  likely  lies  somewhere  between  these  two 
limits. 

In  Fig.  2.1(a),  the  cross-hatching  represents  material  in  the 
atmosphere  of  the  primary.  The  atmosphere  has  expanded  to  completely 
fill  the  critical  potential  lobe.  Under  these  conditions,  material 
can  now  flow  onto  the  secondary,  since  it  is  at  the  same  potential  as 
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the  primary.  Angular  momentum  causes  the  material  to  spiral  in 
around  the  secondary,  creating  an  accretion  disk.  Material 
transferring  from  the  inner  edge  of  the  disk  to  the  surface  of  the 
secondary  releases  gravitational  potential  energy,  thus  powering  the 
X-ray  luminosity.  A  shocked  region  may  also  contribute  to  the 
radiation  at  the  point  where  the  inflowing  stream  joins  the  outer 
edge  of  the  accretion  disk. 

Evolutionary  senarios  developed  to  explain  the  formation  of 
massive  binary  stars  (not  necessarily  X-ray  binaries)  have  shown  that 
the  necessary  mass-exchange  rates  are  quite  easily  obtained  (Kopal 
1959;  Paczynski  1971).  Potential  overflow,  however,  presents  us  with 
the  major  problem  of  too  rapid  an  exchange  of  the  material.  The 
exchange  rate  tends  to  speed  up  as  the  primary  loses  mass  and  the 
secondary  gains  it.  In  a  very  short  period  of  time,  the  material 
flow  can  completely  smother  the  secondary  in  a  X-ray  opaque  cloud. 

The  second  mechanism  is  gravitational  capture  of  material  from 
the  primary's  stellar  wind.  Hoyle  and  Lyttleton  (1939)  developed  the 
analytical  line  accretion  model  to  explain  how  a  galaxy  moving 
through  the  intergalactic  medium  could  gain  matter.  They  argued  that 
the  mass  accretion  rate  is  dependent  only  on  the  density  of  the 
passing  material,  its  velocity,  and  the  mass  of  the  gravitating  body. 
An  accretion  radius  is  defined  by 

R  =  2  G  M  v"2  ,  (2.4) 

a 

where  G  is  the  gravitational  constant,  M  is  the  mass  of  the 
gravitating  body,  and  v  is  the  velocity  the  material  would  have  if 
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the  gravitating  body  were  not  present.  The  unperturbed  kinetic 
energy  of  material  within  this  radius  is  less  than  the  gravitational 
potential  energy.  Under  such  conditions,  we  have  a  simple  criterion 
for  determining  whether  or  not  the  material  can  be  gravitationally 
captured.  Since  material  within  this  radius  is  energetically  favored 
to  be  captured,  the  capture  rate  is  dependent  only  on  the  mass  flux 
through  the  cross  section  defined  by  R  and  some  efficiency  factor. 

a 

So,  we  find  the  accretion  rate  to  be  given  by 

dM/dt  =  it  P  v  a  ,  (2.5) 

a 

where  p  is  the  unperturbed  density  of  the  wind  near  the  body,  and  <* 

is  an  efficiency  parameter  lying  between  0.5  and  1.  Hunt  (1971, 

1979)  performed  numerical  calculations  of  gravitating  bodies  moving 

through  an  interstellar  gas  at  various  Mach  flow  numbers.  He 

demonstrated  that  a  tends  to  increase  with  increasing  Mach  number. 

This  increase  is  attributable  to  the  fact  that  proportionately  more 

energy  is  dissipated  in  the  shock  at  high  Mach  numbers  than  at  lower 

ones.  Thus  material  which  has  passed  through  a  highly  supersonic 

shock  is  easier  to  capture  gravitationally.  Davidson  and  Ostriker 

combined  this  theory  with  the  observations  that  0  and  B  class  stars 

-7  -4  -1 

have  high  mass-loss  rates  on  the  order  of  10  to  10  M0  yr 
(Weymann  1963;  Morton  1967a, b,  1976;  Conti  and  Cowley  1975;  Hutchings 
1976;  Conti  1978a,  Conti  and  Germany  1980).  They  then  noted  that 
these  mass  loss  rates  can  provide  winds  of  suitable  density  and 
velocity  near  the  secondary  so  as  to  allow  the  line-accretion  model 
to  support  the  observed  luminosity.  They  therefore  concluded  that 
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gravitational  accretion  from  a  stellar  wind  could  be  one  workable 
method  of  powering  the  binary  X-ray  stars.  This  mechanism  is 
depicted  in  Fig.  2.1(b). 

Gravitational  accretion  is  also  the  most  likely  method  for 
powering  the  possible  third  class  of  X-ray  binary.  In  these  systems, 
the  secondary  is  pictured  as  orbiting  within  the  extended  tenuous 
atmosphere  of  the  primary  (Amnuel,  Guseinov,  and  Rakhamimov  1979  and 
sources  therein).  At  the  present  time  there  is  insufficient 
observational  data  on  these  systems  to  determine  the  actual  powering 
mechanism.  Binary  system  parameters  are  needed  to  determine  the 
velocities  of  the  secondaries  and  the  density  of  the  atmospheres  in 
their  orbits. 

In  either  of  the  two  methods  discussed  above,  there  is  a  limit 
to  the  maximum  attainable  X-ray  luminosity.  This  is  the  so  called 
Eddington  luminosity,  or  Eddington  limit,  at  which  the  radiation 
pressure  due  to  the  emergent  X-rays  balances  the  force  of  gravity. 
In  the  case  of  spherically  symmetric  accretion  in  which  Thompson 
scattering  is  the  dominant  opacity,  this  limiting  luminosity  for 
material  of  cosmic  abundances  is  given  by  (Eddington  1926) 


'edd 


=  4  u  G  M  c 


(2.6) 


or, 

Ledd  “  1*3xl°38  M/M0  er9s  sec'1,  (2.7) 
where  kq  is  the  electron  scattering  opacity.  Any  excess  material  in 
the  accretion  flow  is  blown  off  by  the  radiation  pressure.  Detailed 
one-dimensional  models  of  spherical  accretion  near  the  Eddington 
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limit  have  been  presented  by  Vi  tel lo  (1978).  He  verified  that  sharp 
cut-offs  near  the  Eddington  limit  existed  for  several  types  of  flow. 
He  also  speculated  that  any  instabilities  introduced  into  the  flow 
may  result  in  luminosities  higher  than  the  Eddington  limit.  Two 
dimensional  instabilities,  such  as  the  Rayleigh-Taylor  instability, 
or  the  creation  of  photon  bubbles  as  suggested  by  Prendergast  and 
Spiegel  (1973),  may  provide  the  necessary  instabilities  in  the  flow. 

C.  System  Life  Time  Considerations 

We  are  now  faced  with  the  problem  that  both  of  these  mechanisms 
are  restricted  in  their  applicability.  The  potential  overflow 
mechanism  appears  to  be  too  effective  at  providing  material,  while 
accretion  from  a  stellar  wind  does  not  work  well  for  primaries  with 
masses  less  than  19  M0.  Van  den  Heuval  (1975)  demonstrated  the 
limits  of  applicability  for  each  of  these  mechanisms  by  considering 
the  lifetimes  of  various  binary  systems  as  a  function  of  their 
masses.  His  argument  proceeds  as  follows. 

The  time  scale  over  which  the  primary  will  lose  80%  of  its  mass 
is  given  by 

t  .  3xl07  (H0/Mg)2  Rg/R0  Lg/L0  ,  (2.8) 
with  t  given  in  years.  From  this  mass  loss  rate.  Van  den  Heuval 
arrived  at  an  average  accretion  rate  given  by 

dMx/dt  =  2.66xl0~8  Ro/R0  LQ/L0  MQ/M0  ,  (2.9) 
with  the  accretion  rate  given  in  M0  yr“*.  If  the  accretion  rate  onto 
the  secondary  exceeds  about  10“8  M0  yr’*,  then  the  secondary  will 
become  engulfed  in  an  optically  thick  cloud,  allowing  no  X-rays  to 
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escape.  However,  there  is  a  period  at  the  onset  of  the  overflow 
during  which  the  mass  loss  rate  and  the  consequent  accretion  rate  of 
Eq.  2.9  is  much  lower.  This  period  exists  for  approximately  1%  of 
the  total  exchange  time  given  by  Eq.  2.8. 

From  Eq.  2.8  and  2.9,  Van  den  Hueval  argued  that  systems  with  a 
primary  more  massive  than  2.1  Mg  would  produce  an  overflow  rate 
sufficient  to  smoother  the  X-ray  source.  In  addition,  the  brief 
start-up  period  from  Eq.  2.8  would  be  on  the  order  of  1000  years  or 
less.  This  brief  period  would  make  it  highly  unlikely  that  such  a 
source  could  be  observed  (see  also  Savonije  1979;  Thomas  1977). 

The  case  for  systems  with  primaries  less  massive  than  2.1  Mg  is 
different.  The  peak  accretion  rates  given  by  Eq.  2.9  are  less  than 
10“6  Mg  yr“*.  Thus  the  X-ray  radiation  is  not  likely  to  be 
smoothered,  and  the  mass-exchange  period  is  very  long.  Hence,  it's 
much  more  probable  that  such  a  system  could  be  observed. 

As  for  accretion  from  a  stellar  wind,  van  den  Heuvel  argued  that 
primaries  less  massive  than  20  Mg  would  have  wind  densities  too  low 
to  supply  the  needed  accretion  material.  However,  those  more  massive 
than  20  Mg  would  have  sufficient  wind  density.  In  addition,  such 
stars  could  be  expected  to  sustain  these  winds  well  in  excess  of  106 
years.  This  "life  time"  criterion  for  X-ray  binaries  appears  to 
agree  well  with  the  binaries  for  which  we  have  sufficient  system 
data.  So  far,  all  X-ray  binaries  meet  one  of  these  two  criterea. 
Thus,  one  of  the  constraints  on  designing  a  model  driven  by  accretion 
from  a  stellar  wind  is  to  require  the  primary  to  have  a  mass  of  20  Mg 
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or  greater.  This  is  in  addition  to  the  goal  of  attempting  to  model 
the  X-ray  production  so  as  to  be  in  reasonable  agreement  with  the 
observations. 

There  are  recent  optical  observations  of  the  massive  X-ray 
binaries  (Hutchings  1976,  Conti  1978b)  which  appear  to  indicate  that 
the  primaries  fill  or  nearly  fill  their  potential  lobes.  This  fact 
could  be  troublesome  if  the  primaries  do  indeed  fill  the  lobes,  for 
as  noted,  there  is  no  way  to  prevent  the  onset  of  run-away  overflow. 
If  the  primaries  only  approach  filling  their  lobes,  there  may  be  a 
much  more  complicated  accretion  process  going  on.  Alme  and  Wilson 
(1976)  have  shown  that  X-ray  heating  of  the  atmosphere  of  a  primary 
which  is  within  80%  of  filling  its  lobe  can  cause  an  outflow  of 
material  sufficient  to  power  the  X-ray  source.  The  massive  binaries 
may  therefore  be  utilizing  both  accretion  from  the  primaries'  massive 
winds  and  X-ray  induced  overflow.  If  the  wind  in  a  given  system  is 
sufficient  for  powering  the  X-ray  source,  there  may  then  be  a  new 
problem  of  limiting  the  total  accretion  rates  obtained  from  several 
exchange  processes  operating  simultaneously.  Such  circumstances  are 
not  examined  in  this  study. 
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III.  THE  MODELING  METHOD 

One  of  the  major  difficulties  encountered  in  modeling 

hydrodynamic  and  radiation  effects  is  the  large  difference  in  the 

time  scales  over  which  the  two  processes  occur.  The  hydrodynamic 

time  scale  is  governed  by  the  typical  scale  lengths  of  the  problem 

and  the  flow  velocities  (or  the  sound  velocities  if  the  flow  is 

subsonic).  In  the  region  of  the  model  near  the  secondary,  a  typical 

g 

cell  dimension  is  on  the  order  of  10  cm.  Here  the  flow  reaches  its 
highest  velocities  of  about  10  cm  sec'  .  This  requires  that  the 
model  be  able  to  resolve  events  in  the  fluid  flow  which  occur  on  a 
time  scale  of  about  10  seconds.  Conversely,  the  X-ray  heating  times 
are  dominated  by  the  Compton  heating  time  for  the  high  temperatures 
expected  near  the  secondary.  Since  the  Compton  heating  time  will  be 
discussed  more  fully  in  Section  C  of  this  Chapter,  we  will  simply 
note  that  from  Eq.  3.26,  we  have  a  typical  Compton  heating  time  on 
the  order  of  10  sec.  To  explicitly  follow  the  X-ray  heating,  we 

_3 

would  require  time  resolution  of  less  than  10  seconds.  Since 
material  requires  about  25000  seconds  to  move  through  the  problem 
mesh,  such  a  small  time  resolution  would  be  prohibitively  expensive 
in  terms  of  computational  resources.  In  the  developments  which 
follow,  we  describe  a  method  by  which  the  model  can  utilize  a  time 
resolution  suitable  for  the  fluid  flow  and  still  approximate  the 
X-ray  heating  effects. 

A.  The  Equations  of  Motion 

We  begin  our  discussion  of  the  development  of  our  model  systems 
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with  a  presentation  of  the  fluid  equations.  These  were  taken  in 
conservative  form  (Potter  1973)  with  two  additional  terms 
incorporated  to  account  for  X-ray  heating  and  the  stellar  wind 
driving  force.  For  the  conservation  of  mass,  we  have 

3p/3t  +  V-pv  =  0  ,  (3.1) 

and  for  the  conservation  of  energy 

3E/3t  +  V-Ev  =  -P(Vv)  +  Q(V-v)  +  (E$$  -  Ejt"1  ,  (3.2) 

while  the  conservation  of  momentum  is  given  by 

3jD/3t  +  7*£v  =  ptf$  +  pV<f  -  VP  -  7Q  .  (3.3) 

9  w 


where 


p  =  mass  density, 

£  =  momentum  density, 

E  =  energy  density, 

P  =  pressure, 
v,  *  velocity, 

Q  =  Richtmyer-von  Neumann  artifical  viscosity, 

=  gravitational  potential, 

$w  =  effective  wind  potential, 
t  =  effective  X-ray  heating  time, 

Ess  =  steady-state  ener9Y  tfie  gas. 

The  Richtmyer-von  Neumann  artifical  viscosity  is  a  numerical  aid 
required  by  the  finite  difference  equations  used  in  the  computer 
simulation  (see  Richtmyer  and  Morton  1967,  or  Potter  1973).  It 
serves  to  smooth  out  shocks  which  form  in  the  flow  over  several  zones 
of  the  model  mesh,  thus  allowing  the  code  to  approximate  the 
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discontinuity  presented  by  the  shock.  This  will  be  discussed  in 
greater  detail,  along  with  the  differencing  method,  in  Section  D. 

B.  The  Wind  Force  Term 

Models  aimed  at  explaining  the  huge  mass  losses  observed  in  0 
and  B  class  supergiants  have  been  presented  by  Lucy  and  Solomon 
(1970)  and  Castor,  Abbott,  and  Klein  (1975).  Castor,  Abbott  and 
Klein  demonstrated  that  the  source  of  the  forces  driving  the  mass 
loss  was  line-absorption  in  CNO  series  elements,  primarily  the  UV 
lines  produced  by  the  L-shell  electrons.  The  photons  absorbed  by  the 
atoms  are  primarily  traveling  radially  outward  from  the  star.  The 
excited  atoms  then  re-emit  photons  isotropically  as  they  return  to 
their  unexcited  state.  This  results  in  a  net  outward  momentum  gain 
for  the  atoms. 

Castor,  Abbott  and  Klein  took  the  force  to  be  expressed  by 

f  =  M(t)  F  ae  c’1,  (3.4) 

where  ag  is  the  electron  scattering  coefficient,  F  is  the  total 
radiant  flux,  and  M(t)  is  a  force  multiplier  dependent  on  atomic 
populations,  oscillator  strengths,  and  optical  depth.  They  then  set 
up  the  fluid  equations,  assuming  spherically  symmetric  outflow,  and 
specifically  accounted  for  gravity,  gas  pressure,  continuum  radiation 
pressure,  and  line  radiation  pressure.  They  found  a  solution  which 
yields  a  nearly  constant  ratio  of  wind  acceleration  to  gravitational 
acceleration.  The  forces  were  also  found  to  be  sufficient  to  support 
both  the  high  mass  loss  rates  and  the  terminal  wind  velocities 
observed  in  supergiants.  In  particular,  they  were  able  to  determine 
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an  approximate  relationship  for  the  velocity  of  the  wind  as  a 
function  of  its  distance  from  the  primary.  This  velocity  law  is 
given  by 

v  =  v,  (1  -  R0/r)%  ,  (3.5) 
where  vro  is  the  wind  velocity  at  an  infinite  distance  from  the 
primary,  Rq  is  the  wind  initiation  radius  on  the  order  of  the  primary 
radius,  and  r  is  the  radius  of  interest.  The  success  of  their  model 
indicates  that  the  gravitational  accretion  driven  binary  X-ray  source 
model  contains  a  serious  complication.  McCray  (1974)  and  McCray  and 
Hatchett  (1975)  speculated  that  if  the  X-ray  luminosity  was  high 
enough,  the  atoms  responsible  for  the  wind  force  would  be  fully 
stripped  of  their  L-shell  electrons,  suppressing  the  wind  in  those 
regions  so  illuminated.  They  believed  that  if  the  velocity  was 
reduced,  accretion  may  be  facilitated,  but  also  indicated  that  the 
exact  effect  was  uncertain. 

We  note  here  that  they  also  claimed  that  the  orbital  periods 
were  comparable  with  the  time  it  takes  the  wind  to  move  from  the 
primary  surface  to  the  secondary's  orbit,  indicating  that  there  would 
be  strong  azimuthal  asymmetry  in  the  wind  with  respect  to  the 
secondary's  position.  A  simplified  ballistic  analysis  presented  in 
Chapter  VII  found  that  the  time  it  takes  the  wind  to  move  from  the 
surface  of  the  primary  to  the  orbit  of  the  secondary  is  actually  on 
the  order  of  10%  of  the  orbital  period.  Our  Model  3,  which  is 
presented  in  Chapter  VIII,  directly  addresses  the  question  of  what 
impact  the  X-rays  have  on  the  structure  of  the  wind. 
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To  account  for  the  stellar  wind  and  its  driving  force,  we  use 
Eq.  3.5  to  derive  a  potential  for  the  driving  force.  We  first  apply 
the  chain  rule  to  Eq.  3.5,  yielding 

3v/3t  =  3v/3r  3r/3t  .  (3.6) 
Since  3r/3t  is  just  v,  we  can  differentiate  Eq.  3.5  with  respect  to  r 
and  substitute  the  result  into  Eq.  3.6,  giving  us 

3v/3t  =  h  v2  R  /  r2  .  (3.7) 

co  0 

To  get  the  potential,  we  now  set  Eq.  3.7  equal  to  the  negative 
gradient  of  the  potential  so  that  integrating  with  respect  to  r  gives 
the  desired 

v  •  (3-8) 

This  is  in  effect  a  "negative  mass"  gravitational  potential  and  can 
thus  be  directly  incorporated  into  the  conservation  of  momentum 
equation  (Eq.  3.3). 

C.  The  X-ray  Heating  Term 

Of  more  importance  is  the  last  term  of  Eq.  3.2,  representing  the 
X-ray  heating  effect.  Several  previous  works  (Tarter,  Tucker,  and 
Salpeter  1969;  Hatchett,  Buff,  and  McCray  1976)  have  shown  that  the 
ionization  and  temperature  structure  of  an  optically  thin  gas  with 
cosmic  abundances  can  be  completely  described  in  terms  of  a  single 
parameter,  £,  for  a  given  spectral  shape,  where 

E  =  Lx  /  n  r2  .  (3.9) 
Here,  Lx  is  the  X-ray  luminosity,  n  the  number  density  and  r^  the 
distance  from  the  X-ray  source. 

Using  e,  as  the  independent  parameter,  Tarter  (see  reference  in 
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Alme  and  Wilson  (1974))  and  Alme  and  Wilson  (1974)  calculated  the 
resultant  equilibrium  temperatures  for  an  optically  thin  gas  exposed 
to  X-rays.  Since  the  results  are  dependent  on  the  spectrum  of  the 
X-rays,  they  both  used  a  22-keV  exponential  spectrum  representative 
of  cosmic  X-ray  sources.  The  results  of  their  calculations  are 
displayed  as  curves  1  and  2  in  Fig.  3.1.  Figure  3.1  thus  allows  the 
equilibrium  temperature  of  a  gas  to  be  specified  for  any  given  £. 
Curves  3  and  4  serve  as  analytical  fits  for  just  this  purpose  and 
will  be  utilized  later  in  the  determination  of  the  X-ray  heating 
contribution  to  the  hydrodynamic  equations. 

For  compatibility  with  the  computer  code  adapted  for  the  models, 
we  have  defined  K  to  be 

5  =  P  r2x  /  Lx  .  (3.10) 
We  can  relate  Eqs.  3.9  and  3.10  through  p  and  n  by  first  noting  that 
(Clayton  1968) 

P  =  ny/NQ  ,  (3.11) 
where  Nq  is  Avogadro's  number  and  y  is  the  mean  atomic  weight.  Now, 
y  is  given  by 

V  =  l  fznz/Az  ,  (3.12) 
where  fz  is  the  fraction  by  weight  of  atomic  element  Z,  nz  is  the 
number  of  particals  contributed  by  element  Z,  and  Az  is  the  atomic 
weight  of  element  Z.  For  an  un-ionized  gas  of  cosmic  abundances,  we 
have 

y  =  l  fz/Az  =  1.16198  .  (3.13) 
Thus,  we  find  that  p  and  n  are  related  through  Eq.  3.12  by 
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Figure  3.1.  Steady  State  Temperature  vs  £.  Steady  state  temperature 
as  a  function  of  t,  is  displayed  here  for  an  optically  thin  gas 
exposed  to  a  22  keV  exponential  X-ray  source.  Curve  1  was  produced 
by  Tarter  (see  Alme  and  Wilson  1974),  while  Curve  2  was  produced  by 
Alme  and  Wilson  (1974).  Curve  3  is  an  analytical  fit  to  the  upper 
portion  and  is  given  by  Eq.  19.  Curve  4  is  a  fit  to  the  tail  and  is 
given  by  Eq.  20.  C  is  not  a  dimensionless  constant,  but  rather  has 
units  of  sec  3  cm'3. 


P  =  1.9294xl0‘24  n. 


(3.14) 


The  X-ray  heating  rate,  including  photoionization  and  Compton 


scattering,  is  given  by: 


(3E/3t)r  =  (E$s  -  E)  t;' 


(3.15) 


Here,  t  is  the  X-ray  heating  time,  and  E  is  the  steady  state 

A  bb 

(equilibrium)  energy  of  the  gas.  E$$  is  computed  using 

Ess  =  p  Cv  Tss  *  (3-16) 

with  Tss  taken  from  Fig.  3.1  for  a  given  K. 

Two  analytical  fits  were  developed  to  give  the  steady  state 

temperature  as  a  function  of  K.  The  upper  portion  of  the  curve  is 

fitted  with  the  dashed  line  (curve  3)  specified  by  (T  in  ev) 


1  +  7.716xl054  K2 


(3.17) 


The  roll-off  region  of  the  tail  is  fitted  with  the  straight  line 


(curve  4)  specified  by  (T  in  ev) 


T  =  1.5689x10’ 9  C0-362704 


(3.18) 


From  Eq.  3.15,  we  see  that  if  the  gas  is  below  its  steady  state 
temperature,  heating  will  occur,  while  if  it  is  above  its  steady 
state  temperature,  then  cooling  will  occur.  The  time  over  which  this 
can  occur  is  called  the  X-ray  heating  time.  We  can  estimate  the 
X-ray  heating  time  by  first  calculating  the  Compton  heating  time. 


This  time  can  be  given  by 

tc  =  E  (dE/dt)'1 
where  ( dE/dt ) c  is  the  Compton  heating  rate. 


(3.19) 


The  Compton  heating  rate  in  an  optically  thin  gas  can  be  derived 
from  one  term  of  the  radiation  diffusion  equation  in  the  form  (Alme 


and  Wilson  1974) 

(3E/3t)c  =  p/mec  <e(kT  (v  3Ev/3v  -  3Ey)  +  hvEv>  dv  ,  (3.20) 
where  induced  scattering  has  been  neglected.  We  have  assumed  an 
exponential  spectrum  for  our  models,  so  we  have 

■  A  •  (3-21) 
where  A  is  chosen  so  as  to  give  the  magnitude  of  the  spectrum  and  vq 
is  the  characteristic  frequency  of  the  spectrum.  In  this  case,  hvQ 
is  22  keV. 

Substituting  Eq.  3.21  into  3.20  and  integrating  yields 

(3E/3t)c  =  (hvQ  -  4kT)  P  <e  vq  A  /  me  i.  .  (3.22) 
Since  the  integral  of  Eq.  3.21  over  all  v  gives  the  total  energy  in 
the  spectrum,  we  must  have  vqA  equal  to  the  radiation  energy  density, 
E  .  Equation  3.22  specifies  the  conditions  under  which  the  material 
will  come  into  equilibrium  with  the  radiation.  We  are  now  able  to 
define  a  characteristic  Compton  heating  time  by 

(3E/3t)c  =  Er  P  <e  <hv>  /  c  me  ,  (3.23) 
where  <hv>  is  now  the  characteristic  equilibrium  energy  due  to  the 
X-ray  spectrum. 

The  internal  energy  of  the  gas  can  be  given  by 

E  =  P  Cy  T  .  (3.24) 
We  also  have  the  radiation  energy  density  due  to  the  X-ray  luminosity 
given  by 

Er  =  Lx  '  4  11  r2  c  *  (3.25) 
Dividing  Eq.  3.24  by  Eq.  3.23  and  substituting  Eq.  3.25  for  Er,  we 
arrive  at  a  characteristic  Compton  heating  time  given  by 
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(3.26) 


4  tt  r2  m  C ,  T 

t  =  - §_ ~ — 

c  L  <  <hv> 
x  e 

For  gas  of  low  temperature,  the  photoelectric  effect  dominates 
the  X-ray's  interaction  with  the  gas.  The  ionization  of  atoms  in  the 
gas  effectively  absorbs  energy  from  the  incident  radiation.  Above 
some  transition  energy  E^,  Compton  scattering  becomes  the  dominate 
energy  transfer  process  between  the  radiation  and  the  gas.  The  gas 
has  become  highly  ionized  and  the  atoms  are  no  longer  able  to  absorb 
energy  by  the  photoelectric  effect.  The  X-ray  heating  time  is 
therefore  approximated  by 


t 


x 


1  +  (Et/E)2  ' 


(3.27) 


Equation  3.27  allows  the  X-ray  heating  time  to  approach  the  Compton 

time  whenever  E  is  greater  that  Efc,  and  strongly  magnifies  the  X-ray 

heating  effects  for  E  less  than  E^.. 

The  material  in  the  stellar  wind  is  assumed  to  be  optically 

thin.  This  is  based  on  the  fact  that  the  wind  density  in  the  region 

-12  -3 

of  the  secondary  is  found  to  be  on  the  order  of  10  g  cm  for  mass 
loss  rates  required  to  support  the  observed  X-ray  luminosities. 
Further,  with  the  aid  of  Fig.  3.1  and  Eq.  3.10,  we  note  that  strong 
X-ray  heating  will  occur  only  close  to  the  secondary.  Thus  our  X-ray 
heating  model  will  retain  its  validity  even  if  the  material  becomes 
optically  thick  far  from  the  region  of  strong  heating.  Optical 
depths  were  calculated  through  various  regions  of  the  models  after 
they  were  run  to  determine  the  validity  of  this  assumption.  The 
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resultant  optical  depths,  to  be  discussed  in  Chapter  IX,  show  that 
the  models  are  in  fact  optically  thin  near  the  secondary,  but  became 
progressively  thicker  in  the  outer  regions. 


D.  The  Numerical  Method 

The  equations  of  motion  (Eqs.  3.1-3)  are  solved  utilizing  the 
Lax-Wendroff  two-step  method  for  finite-difference  equations 
{Richtmyer  and  Morton  1967;  Potter  1973).  This  method  defines  the 
dependent  variables  on  alternate  mesh  points  from  the  corresponding 
fluxes.  With  respect  to  Fig.  3.2,  this  can  be  visualized  by  defining 
the  dependent  variables  at  the  center  of  each  cell  and  the  fluxes  at 
the  center  of  each  face  of  the  cell.  The  two-step  method  first 
computes  the  new  fluxes  from  the  old  dependent  variables: 

Uj+%  =  %  (uj  +  Uj+1}  '  (Fj+l  "  Fj)  At  /  2  Ax  .  (3.28) 

This  is  then  used  to  give  the  fluxes  from 

(3.29) 

Here,  n  refers  to  the  n-th  time  step  while  j  refers  to  the  j-th 
spacial  increment.  The  new  dependent  variables  are  then  given  by 
these  fluxes  as 


Pn+Js  =  F(un+?S) 

hj-ns  lj+v 


u3+1  =  u3  ■  (FjS "  FjS}  At/Ax  •  (3-30) 

We  are  now  in  a  better  position  to  discuss  the  requirement  for 
the  Richtmyer-von  Neumann  artificial  viscosity,  Q,  introduced  in  Eqs. 
3.2  and  3.3.  From  the  finite  dimensions  utilized  by  the  differencing 
scheme  described  above,  it  is  apparent  that  only  phenomena  with 
wavelengths  greater  than  the  dimensions  of  the  cells  will  be 
accurately  described.  In  a  compressible  fluid,  large  amplitude 


k-l,j-l 


k.j-1 


Figure  3.2.  The  Lax-Wendroff  Mesh.  Shown  here  is  the  mesh  used  for 
the  Lax-Wendroff  two-step  finite  differencing  scheme.  The  dependent 
variables  p  and  E  are  defined  on  the  "main"  mesh  (dashed  lines), 
while  the  momenta  and  other  fluxes  are  defined  on  the  auxiliary  mesh 
(solid  lines).  The  two  meshes  are  most  easily  pictured  as  one  set  of 
cells  with  variables  defined  as  zone-centered  or  face-centered. 
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resolution  near  the  secondary  and  a  large  mesh  size  with  a  limited 

g 

number  of  zones.  A  minimum  cell  dimension  of  2.5x10  cm  represents  a 
compromise  between  resolution  needed  near  the  accretion  region  of  the 
secondary  and  computational  time  available.  Since  the  primary  thrust 
of  this  study  is  to  examine  the  flow  characteristics  and  not  just  the 
details  of  accretion  we  do  not  feel  that  this  lack  of  detail  close  to 
the  secondary  is  significant. 

While  the  wind  flow  is  more  correctly  described  as  a  plasma 

rather  than  as  a  gas,  magnetic  effects  have  been  omitted  in  these 

models.  One  reason  is  that  the  minimum  dimension  of  our  models  is 

7  8 

well  outside  the  10  to  10  cm  magnetospheric  radii  of  neutron  stars 
(Lamb,  Fabian  and  Pringle  1973).  Another  is  that  there  is  little 
information  as  to  the  effects  of  OB  supergiant  magnetic  fields  on 
their  massive  winds.  In  addition,  Parker  (1972)  has  presented  a 
mechanism  by  which  such  fields,  if  entrained  in  the  wind,  could  be 
dissipated.  The  final  reason  is  that  treatment  of  the  magnetic 
fields  would  increase  the  difficulty  of  the  problem  to  the  point 
where  it  could  not  be  handled  with  the  computational  resources 
available. 

We  have  limited  the  models  to  two  dimensions  for  several  reason. 
The  first  is  a  practical  one,  namely  the  lack  of  the  computational 
resources  necessary  for  a  full  three  dimensional  model  of  this  type. 
We  also  feel,  however,  that  it  is  desirable  to  maintain  a  link 
between  previous  two-dimensional  models  that  did  not  include 
radiation  effects  and  the  present  one.  A  second  consideration  is 
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that 

the 

time-scales 

for  the 

flows  near  the 

accretion  region  are  an 

order 

of 

magnitude 

faster 

than  the 

orbital 

motion  time  scales,  so 

that 

the 

exclusion 

of  the 

orbital 

motion 

does  not  significantly 

impact  on  the  accretion  flow.  The  lack  of  orbital  motion  does 
prevent  any  possible  formation  of  an  accretion  disk  around  the 
secondary  since  the  angular  momentum  of  the  wind  is  ignored.  In 
Chapter  VII,  we  discuss  some  of  the  consequences  brought  on  by  the 
three-dimensional  nature  of  the  actual  flows. 
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IV.  GENERATING  MODEL  PARAMETERS 

In  the  previous  chapter,  we  outlined  the  modeling  method  to  be 
used.  The  next  step  is  to  establish  the  parameters  for  our  model 
binary  systems.  The  parameters  are  selected  to  enforce 
self-consistency  between  the  accretion  rate  and  the  desired  model 
luminosity.  To  achieve  this  goal,  we  need  to  determine  how  the  wind 
density  and  velocity  vary  at  the  orbit  of  the  secondary  with  changes 
in  the  other  system  parameters. 

We  begin  our  analysis  by  deriving  a  relationship  for  the 
luminosity  as  a  function  of  the  stellar  wind  density  and  velocity. 
Let  us  again  take  the  wind  velocity  to  be  given  by  Eq.  3.5  as 

v  =  v.  (1  -  \/r)'2  .  (4.1) 

We  also  adopt  the  results  of  Zel'dovich  and  Shakura  (1969),  Shakura 
and  Sunyaev  (1973),  and  Alme  and  Wilson  (1973)  and  assume  that  the 
gravitational  potential  energy  released  by  the  infalling  matter  is 
equivalent  to  102!  of  the  matter's  rest  mass  energy.  The  luminosity 
in  terms  of  the  accretion  rate  is  then  given  by 

Lx  =  0.1  dMx/dt  c2  .  (4.2) 

Since  the  secondary's  accretion  rate  is  given  from  Eq.  2.5  as 

dMx/dt  =  it  R2  Px  vx  a  ,  (4.3) 

we  find  the  luminosity  to  be 

Lx  =  0,1  11  Ra  px  vx  a  °2  ’  (4-4) 

where  Px  and  vx  are,  respectively,  the  unperturbed  wind  density  and 

velocity  at  the  orbit  of  the  secondary.  We  can  now  ;se  Eq.  2.4  to 

specify  the  accretion  radius  of  the  secondary: 
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K  =  2  G  v‘2  .  (4.5) 

Substituting  Eq.  4.5  into  Eq.  4.4  we  arrive  at  the  desired  equation 
giving  the  luminosity  as  a  function  of  the  wind  density,  wind 
velocity,  accretion  efficiency,  and  the  mass  of  the  secondary: 

Lx  =  0.4  tt  ce  px  (cGMx)2  v'3  .  (4.6) 

Equation  4.6  points  out  the  strong  dependence  of  the  X-ray  luminosity 
on  the  mass  of  the  secondary  and  the  velocity  of  the  wind. 

The  accretion  rates  found  in  the  models  are  compared  to  those 
predicted  by  Eq.  4.3.  This  defines  a  ratio  between  the  predicted 
accretion  rate  and  the  resulting  rate.  The  accretion  efficiency  in 
Eq.  4.3  is  initially  estimated  from  Hunt's  (1971)  work. 

We  can  estimate  the  variation  in  the  wind  density  and  velocity 
with  the  binary  system  parameters.  To  do  this,  we  assume  that  the 
wind  is  unperturbed  as  it  passes  the  secondary  and  that  the  X-ray 
luminosity  is  held  constant.  We  also  assume  that  the  primary's  mass 
loss  rate  is  constant.  This  gives  the  density  in  the  wind  as  a 
function  of  the  distance  from  the  primary  as 

dM  /dt 

P= - •  (4.7) 

4  n  r  v 

We  now  substitute  Eq.  4.1  into  Eq.  4.7  and  take  r  to  be  the  binary 
separation.  This  gives  us  the  wind  density  at  the  orbit  of  the 
secondary  as  a  function  of  the  primary's  mass  loss  rate  and  the 
terminal  wind  velocity: 

dM  /dt 

P  = - o-2 - c  •  (4-8) 

X  4TtaSoe(l-R  /ap 
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In  our  later  discussions,  it  will  prove  useful  to  have  an  expression 
for  5  as  a  function  of  distance  from  the  primary  and  secondary.  From 
Eq.  3.10  we  have 


5  *  0  rx  /  Lx  • 

We  can  now  incorporate  Eqs.  4.1  and  4.7  into  Eq.  4.9  to  obtain 


(4.9) 


dVdt  rx 


(4.10) 


*  4  tt  r2  L  v  (1-R  /r  )h  ' 

0  X  oo  v  0  0 

where  rQ  is  the  distance  from  the  primary. 

One  of  the  more  critical  determinations  to  be  made  is  the 
primary  mass  loss  rate  needed  to  support  the  secondary's  required 
accretion  rate.  We  derive  this  from  Eq.  4.8  using  Eqs.  4.3  and  4.5. 
This  yields 


dMQ/dt  = 


dMx/dt  vj  (a-Ro)2 
a  (GMx)2 


(4.11) 


And  lastly  the  wind  initiation  radius  as  a  function  of  the  terminal 
wind  velocity,  the  primary's  mass  loss  rate,  and  the  secondary's 
accretion  rate  are  found  from  Eq.  4.11  as 


R0  «  a  -  (GMx/v^)  (a  dM0/dt  /  dMx/dt))5'2  .  (4.12) 
Equations  4.11  and  4.12  determine  our  model  parameters  for  they 
include  the  major  variables  describing  a  binary  system. 

In  establishing  our  model  parameters,  we  desire  reasonably  close 
agreement  with  the  observed  systems.  However,  our  self-consistency 
constraint  may  require  some  adjustment  of  the  model  parameters  away 
from  their  observed  values.  The  three  system  parameters  with  the 
highest  observational  uncertainty  are  the  most  reasonable  choices. 
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They  are  dM  /dt  (or  equivalently,  L  ) ,  R  .  and  vm.  To  a  lesser 
extent,  Mx  could  be  varied,  but  observations  place  tighter 
constraints  on  it.  RQ  and  both  affect  the  wind  velocity  through 
Eq.  4.1,  and  the  sensitivity  of  the  X-ray  luminosity  to  the  wind 
velocity  is  plainly  evident  from  Eq.  4.6. 

In  addition  to  observational  constrains,  there  are  theoretical 
limitations  to  be  considered  in  parameter  selection.  Rq  has  limits 
imposed  on  it  by  the  theory  of  Castor,  Abbott  and  Klein  (1975).  The 
initiation  radius  is  essentially  the  photospheric  radius  of  the 
primary,  although  it  may  be  as  much  as  25%  greater. 

Constraints  are  also  set  on  acceptable  values  of  the  primary's 
mass  loss  rate  from  energy  and  momentum  arguments  (for  review,  see 
Cassinelli  1979).  First,  assume  that  a  photon  is  emitted  in  the 
photosphere  of  the  primary  with  a  frequency  larger  than  some  strong 
line.  If  the  atmosphere  is  expanding  radially  outward  and 
accelerating,  then  it  may  eventually  reach  a  velocity  at  which  the 
photon  is  Doppler  shifted  into  resonance  with  the  line.  This 
velocity  is  given  by 

v  =  c(v-vq)/vo  ,  (4.13) 
where  vq  is  the  line  central  frequency.  For  a  single  line,  we  can 
equate  the  final  mass  momentum  flux,  dMQ/dt,  to  the  photon 
momentum  that  is  transferred  by  scattering  all  the  radiation  between 

vo  and  vo  p1us  V^c-  This  9ives  us 

v.  dM0/dt  =  LvoVa/c2  .  (4.14) 

If  the  entire  spectrum  is  covered  by  non-overlapping  lines  such  that 


44 


each  adjacent  line  is  seperated  by  a  displacement  corresponding  to 
the  Doppler  shift  at  vM,  the  maximum  mass  loss  rate  can  be  related  to 
the  total  momentum  flux  of  the  luminosity  giving 

<dVdt>max  ■  L'*~c  •  <4.15) 
Note  that  the  kinetic  energy  of  the  mass  loss  is  given  by 

H  dMQ/dt  v^  =  *5  L  vffl  c-1  .  (4.16) 
3  -1 

Typical  values  of  v^  are  1-2x10  km  sec  ,  so  the  mass  loss  carries 
away  about  0.5%  of  the  radiative  luminosity  of  the  star.  Castor  (see 
reference  in  Casinelli  1979)  has  shown  that  if  multiple  scattering  of 
the  photons  occurs,  the  mass  loss  given  by  Eq.  4.15  will  be  raised  by 
not  more  than  a  factor  of  two  or  three. 

A  model  is  generated  by  starting  with  the  observed  values  for 
the  parameters  of  the  desired  binary  system.  Using  these  values,  we 
then  calculate  the  required  accretion  rate  using  Eq.  4.2.  Equation 
4.11  is  then  used  to  calculate  the  required  mass  loss  rate  from  the 
primary.  This  is  checked  against  the  maximum  rate  given  by  Eq.  4.15 
to  determine  its  acceptability.  If  the  required  value  is  not 
acceptable,  one  or  more  of  the  parameters  are  adjusted  and  the  above 
process  is  repeated  until  we  arrive  at  an  acceptable  value  for  the 
mass  loss  rate.  At  this  point,  the  resulting  set  of  parameter  values 
gives  us  our  desired  self-consistent  model. 
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V.  MODEL  ONE 

The  primary  purpose  of  Model  1  is  to  provide  a  link  between  the 
two-dimensional  hydrodynamic  studies  of  the  past  and  the 
hydrodynamics  with  radiation  and  wind  ionization  effects  of  Model  3. 
Thus,  Model  1  does  not  account  for  X-ray  ionization  effects  on  the 
wind  force.  It  does,  however,  include  the  X-ray  heating  and  wind 
force  models  discussed  in  Chapter  III.  Modification  of  the  wind 
force  due  to  X-ray  ionization  will  be  left  for  Model  3. 

A.  Initializing  Model  1 

As  a  basis  for  Model  1,  we  select  Vela  X-l  (4U0900-40).  This 
system  is  commonly  assumed  to  be  powered  predominantly  by  wind 
accretion  (Petterson  1978).  It  also  has  one  of  the  largest  X-ray 
binary  separations  known.  The  large  separation  serves  both  as  a 
dominant  factor  in  determining  the  mode  of  powering  the  system  and  as 
an  aid  in  establishing  our  model.  The  zoning  scheme  for  Model  1  does 
not  provide  sufficient  spatial  resolution  to  accurately  represent  the 
strong  wind  acceleration  near  the  primary.  Yet  at  the  same  time,  we 
require  a  large  area  about  the  secondary  in  order  to  minimize 
boundary  effects.  The  Vela  X-l  system  is  best  suited  to  the 
satisfaction  of  these  conditions. 

The  model  parameters  for  this  system  (and  the  4U1700-37  system 
of  Model  3)  are  adjusted  to  be  in  reasonable  agreement  with  the 
observations  reported  by  Avni  and  Bahcall  1975;  Becker,  et  al  1978; 
Greenstein  and  McClintock  1976;  McClintock,  et  al  1976;  Petterson 
1978;  Pravdo,  et  al  1976;  Rappaport,  Joss  and  McClintock  1976; 
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Bahcall  1978;  Hutchings  1974;  and  Ulmer,  et  al  1972  (re.  Table  2.1). 

At  the  same  time,  we  also  seek  to  achieve  the  desired 

self-consistency  between  the  accretion  rate  and  the  luminosity. 

In  Chapter  IV,  we  presented  the  tools  required  to  generate  the 

model  parameters.  For  Model  1,  we  begin  with  the  observed  luminosity 
37  -1 

of  10  erg  sec  .  Since  we  have  the  luminosity  given  in  terms  of 
the  accretion  rate  by  Eq.  4.2,  it  is  a  simple  matter  to  rewrite  it  in 
a  form  giving  the  accretion  rate  in  terms  of  the  desired  luminosity: 

dMx/dt  =  10  Lx/c2  .  (5.1) 

For  the  observed  luminosity,  Eq.  5.1  yields  a  required  accretion  rate 
of  1.1x10*7  g  sec"*. 

In  order  to  estimate  the  accretion  efficiency,  we  note  that  far 
from  the  X-ray  source  the  gas  can  be  expected  to  be  highly 
supersonic,  while  closer  in  it  may  become  subsonic.  For  the  purposes 
of  Model  1,  we  fix  ot  at  0.75,  which  corresponds  to  a  flow  of  Mach  2 
and  is  the  median  value  of  a  between  subsonic  and  highly  supersonic 
flow.  Later  analysis  of  the  flow  will  show  that  this  is  a  reasonable 
value. 


Referring  back  to  Table  2.1,  we  see  that  the  wind  velocity  at 
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infinity  is  1430  km  sec  ,  the  luminosity  of  the  primary  is  1.5x10 

ergs  sec"*,  the  binary  separation  is  52  Rg,  and  the  radius  of  the 

primary  is  33  Rg.  Now,  the  required  primary  mass  loss  rate  is  given 
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With  the  values  cited  this  yields 


dMQ/dt  =  2.49xlO"5  M0  yr_1  .  (5.3) 

To  check  the  feasibility  of  this  value,  we  recall  that  Eq.  4.16 
and  the  discussion  of  Chapter  IV  set  an  upper  limit  on  the  mass  loss 
rate  of 

(dM0/dt)max  .  3  L/v„c  .  (5.4) 

Based  on  the  parameters  for  the  Vela  X-l  primary,  Eq.  5.4  gives  an 
upper  limit  of 

(dMo/dt)max  =  1*66x10'5  M0  ^_1  •  (5.5) 

Hence,  our  required  mass  loss  rate  based  initially  on  the  observed 

parameters  is  larger  than  the  theoretical  maximum  by  a  factor  of 

1.48.  We  elected  to  use  this  value  and  not  pursue  parameter 

modifications  further.  The  model  parameters  are  summarized  in  Table 
5.1.  Note  that  the  mass  loss  rate  is  a  factor  of  1.78  higher  than 
the  upper  limits  set  by  the  observations  of  Hutchings  (1976). 

We  use  these  values  for  several  reasons.  First,  a  48%  increase 
in  the  luminosity  of  the  primary  could  account  for  the  additional 
mass  loss  required.  The  primary  mass  loss  rate  was  also  designated 
as  one  of  the  primary  candidates  for  modification  if  changes  were 
necessary.  Second,  due  to  the  close  initial  match,  we  feel  it  more 
advisable  to  leave  the  majority  of  the  parameters  at  their  accepted 
values  rather  than  modify  them  further.  Third,  we  are  artificially 
forcing  the  luminosity  to  be  sustained  at  its  observed  maximum. 
Model  1  is  not  expected  to  show  any  fluctuation  in  its  accretion 
rates  or  X-ray  luminosity.  Sustaining  the  maximum  luminosity 
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TABLE  5.1.  Parameters  for  Model  1 
Base  System:  4U0900-40  (Vela  X-l) 


Observed 

Model 

a  (Rg) 

52 

52 

w 

33 

33 

v  (km  sec"*) 

1430 

1430 

Mx 

2 

2 

Lx  (ergs  sec"1) 

lxlO37 

lxlO37 

Lq  (ergs  sec"1) 

1.5xl039 

2.2xl039 

dMo/dt  (Mq  yr"1) 

3.5-14.0xl0"6 

2.4886x10 

a 

... 

0.75 

Table  5.1.  l,.e  parameters  used  for  Model  1  are  summarized  here.  The 
observed  values  are  taken  from  Table  2.1.  Note  that  the  only 
difference  in  the  model  parameters  is  the  value  selected  for  the 
primary's  luminosity.  The  selection  of  this  value  is  discussed  in 
the  text. 
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therefore  provides  for  the  largest  possible  X-ray  heating  impact,  and 
avoids  the  selection  of  some  other  arbitrary  value. 

We  initialize  the  model  by  first  establishing  within  the  problem 
mesh  an  unperturbed  wind  flow  with  velocities  given  by  Eq.  4.1  as 

v  =  U  -  >  (5.6) 

and  densities  given  from  Eq.  4.7  as 

dM  /dt 

p  = - o  -  <  (5>7) 

4  11  ro  vx 

The  X-ray  luminosity  is  fixed  at  the  observed  value.  Equations  3.10, 
17,  and  18  are  then  used  to  set  the  material  in  each  cell  to  its 
steady  state  state  temperature. 

We  can  now  examine  the  conditions  within  the  initialized  flow. 
Equation  2.4  for  the  accretion  radius  was 

R  =  2  G  M  v"2  .  (5.8) 

a  X 

Using  the  parameters  for  Model  1,  Eq.  5.8  yields 

R  =  7.09xl010cm.  (5.9) 

a 

At  any  given  point  in  the  wind,  we  can  calculate  the  local  sound 
speed  and  compare  it  with  the  local  wind  velocity.  We  find  that  at 
one  R,  the  flow  is  subsonic  with  a  Mach  number  of  about  0.87.  At  the 

a 

maximum  distance  (8  R,)  from  the  secondary  perpendicular  to  the 

a 

1 ine-of-centers,  we  find  the  flow  to  be  highly  supersonic  with  a  Mach 
number  of  24.  The  flow  first  changes  from  subsonic  to  supersonic  at 
a  distance  of  about  1.4  R  .  This  is  in  agreement  with  our  initial 

a 

arguments  on  the  selection  of  0.75  as  the  accretion  efficiency, 
although  it  is  still  apparent  that  this  value  is  somewhat  arbitrary. 
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The  model  is  then  followed  in  time  until  material  has  flowed 
completely  through  the  spatial  mesh  several  times.  We  note  that  no 
major  changes  occur  in  the  characteristics  of  the  flow  after  this 
time.  The  flow  settles  into  a  quasi-steady  state  such  that  only  very 
slight  fluctuations,  on  the  order  of  0.1%,  are  evident  in  the 
accretion  rate  and  luminosity.  At  this  point  the  calculation  is 
terminated. 

B.  Calculation  Results 

Upon  re-examining  the  conditions  within  the  gas,  we  find  that 
the  supersonic  region  of  the  flow  has  moved  inwards  toward  the 
secondary.  The  flow  remains  highly  supersonic,  nearly  Mach  4.6, 

q 

along  the  line-of-centers  up  to  a  point  5x10  cm  from  the  secondary. 
To  the  sides  of  the  secondary,  the  flow  stays  supersonic  within  R  . 

a 

However,  to  the  rear  of  the  secondary,  in  the  region  of  the  wake,  the 
flow  is  highly  subsonic.  This  behavior  is  related  to  the  temperature 
structure  of  the  gas  as  presented  later  in  this  chapter.  The  flow 
becomes  more  supersonic  than  the  initial  state  indicated.  This 
should  have  resulted  in  a  higher  accretion  ratio  since  the  efficiency 
increases  with  the  Mach  number.  However,  we  find  that  the  accretion 
ratio  actually  decreases. 

In  Chapter  IV,  we  noted  that  Eq.  4.3  can  be  used  to  define  an 
accretion  ratio.  This  ratio,  represented  by  B,  is  given  by 

dM/dt 

B  =  - o-* -  •  (5.10) 

*  Ra  px  vx  a 

The  final  values  for  the  luminosity  and  the  accretion  ratio  reached 


by  Model  1  are 


Lx  =  5.21xl036  erg  sec'1  , 
and 

3  =  0.523  . 

These  results  are  a  factor  of  2  lower  than  those  expected  for 
gravitational  accretion  from  a  planar  flow  as  predicted  by  the  line 
accretion  model. 

We  attribute  the  lower  efficiency  to  the  fact  that  the  wind 
force  is  allowed  to  operate  unaffected  by  the  X-rays.  It  is 
therefore  providing  an  additional  force  acting  to  blow  the  wind  away 
from  the  secondary.  This  force  is  not  included  in  the  line  accretion 
model.  To  see  this,  let  us  consider  the  force  acting  on  the  wind 
along  the  line-of-centers  at  some  distance  r,  which  is  greater  than 
a.  The  expression  for  the  force  is  then 

F  =  Hv3R0/r2  -  GMx/r2  .  (5.11) 
Equation  5.11  goes  to  zero  at  a  distance  of  3.8X1011  cm  from  the 
secondary.  This  is  about  5.4  Rfl.  The  wind  force  is  greater  than  the 
secondary's  gravitational  force  farther  out.  Off  the  line  of 
centers,  the  wind  force  accelerates  the  flow  outwards,  providing  an 
acceleration  component  not  present  in  the  planar  flow  of  the  line 
accretion  model.  Thus,  it  appears  that  the  wind  force  significantly 
affects  the  accretion  flow  in  this  model. 

We  have  produced  a  series  of  figures  which  show  contours  of 
constant  material  density,  energy  density,  temperature  and  5.  There 
are  several  features  in  common  among  these  figures.  The  scales  show 
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distance  from  the  secondary,  with  the  location  of  the  secondary 
marked  by  a  small  x.  In  all  cases,  the  wind  material  is  flowing  in 
from  the  right  and  exiting  to  the  left,  corresponding  to  the  primary 
lying  off  the  right  side  of  the  figure.  The  phenomena  occurring  near 
the  upper  and  lower  edges  of  the  figures  are  boundary  effects.  These 
are  the  result  of  the  spatial  mesh  being  too  small.  The  side 
boundaries  are  not  able  to  supply  material  fast  enough  as  the  flow  is 
drawn  in  by  the  secondary's  gravity.  Hence,  these  regions  become 
artificially  rarified.  This  has  no  effect  on  the  model  since  the 
regions  are  still  well  outside  the  main  accretion  flows. 

Figure  5.1  shows  contours  of  constant  material  density.  The 
phenomenon  occuring  immediately  to  the  right  of  the  secondary  is  the 
result  of  a  thermal  instability  in  the  wind  flow.  To  get  at  the 
nature  of  this  instability,  we  first  note  that  the  energy  contours 
(Fig.  5.2)  are  very  smooth  except  near  the  wake  region  to  the  left  of 
the  secondary.  Since  the  energy  density  is  directly  proportional  to 
the  pressure,  the  smooth  contours  indicate  that  the  temperature  and 
density  are  varying  in  a  manner  which  maintains  a  constant  pressure. 
The  fluctuations  are  wildly  evident  in  Fig.  5.1,  but  are  not  as  clear 
in  the  temperature  contours  of  Fig.  5.3.  This  is  due  to  the  small 
variation  in  density  over  the  mesh,  as  compared  with  a  large 
variation  in  temperature. 

McCray  and  Hatchett  (1975)  demonstrated  that  a  gas  exposed  to 
X-ray  heating  would  exhibit  multiple  density-temperature  states.  A 
gas  under  these  conditions  can  change  states  discontinuously 
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Figure  5.1.  Material  Density  Contours  for  Model  1.  In  this  and  all 
other  contour  plots,  we  have  the  secondary  located  at  the  small  x 
near  the  center  of  the  figure.  The  primary  lies  off  the  right  side 
of  the  plot,  which  is  taken  as  the  "front"  of  the  secondary.  The 
material  flow  is  thus  from  right  to  left  through  the  plot.  In  front 
of  the  secondary  lies  a  region  of  thermal  instability  as  described  in 
the  text.  The  densely  packed  contours  are  at  level  3.  The  contours 
near  the  top  and  bottom  of  the  frame  are  the  result  of  boundary 
effects  in  the  simulation.  The  contour  values  in  10"13  g  cm"3  are: 
(1)  0.5;  (2)  1.0;  (3)  1.6;  (4)  2.6;  (5)  4.1;  (6)  6.6;  (7)  10.0. 
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Figure  5.3.  Temperature  Contours  for  Model  1.  This  plot  clearly 
shows  the  formation  of  a  hot  wake  behind  the  secondary.  Its  regions 
of  instability  are  not  as  evident  in  this  figure  as  in  Fig.  5.1 
because  of  the  much  larger  variation  in  temperatures.  Contours  along 
the  upper  and  lower  edges  are  the  result  of  boundary  effects  in  the 
simulation.  The  contour  values  in  keV  are:  (1)  0.020;  (2)  0.045;  (3) 
0.100;  (4)  0.224;  (5)  0.500;  (6)  1.110;  (7)  2.500. 


depending  on  its  current  state  and  previous  history.  In  particular, 
hot  gases  undergoing  cooling  may  jump  from  a  thin,  high  temperature 
state  to  a  dense,  low  temperature  state,  while  staying  at  constant 
pressure.  Conversely,  a  cool  gas  heating  may  undergo  an  opposite 
transition. 

Following  their  work,  we  have  produced  a  plot  of  P/F  versus  p/F 
for  each  cell  of  the  problem,  where  F,  the  X-ray  flux,  is  given  by 

F  =  - *  ,  (5.12) 

4  tt  r‘ 

and  P,  the  gas  pressure,  is  given  in  terms  of  the  energy  density  of 
the  gas  by 

P  =  E(y-1)  .  (5.13) 

Figure  5.4  clearly  demonstrates  the  behavior  predicted  by  McCray  and 
Hatchett.  Two  different  populations  of  gas  are  evident.  The  upper 
portion  of  the  curve  corresponds  to  hot  gas  which  is  moving  away  from 
the  secondary  and  beginning  to  cool.  The  lower  portion  shows  cool 
gas  which  is  being  heated  as  it  approaches  the  secondary.  The  sharp 
corner  in  the  lower  right  area  of  the  curve  is  an  artifact  of  the 
match  between  Eqs.  3.20  and  3.21  that  were  used  to  fit  the  steady 
state  temperature  curve.  The  gas  attempts  to  change  states 

erratically  in  the  central  region  of  the  plot.  Lines  of  constant 
temperature  drawn  through  the  plot  serve  as  references  for  the 

oc  OC 

multi-valued  region  between  p/F  =  10  and  10  .  The  broken 

contours  in  front  (to  the  right)  of  the  secondary  in  Fig.  5.3 
represent  gas  which  falls  into  this  region  of  the  curve.  As  a  final 


p/F  (sec  cm  ) 


Figure  5.4.  P/F  vs  p/F  for  Model  1.  This  figure  demonstrates  the 
multivaluedness  of  the  pressure-temperature  states  of  the  gas  as  it 
is  heated  or  cooled.  Lines  of  constant  temperature  are  displayed  for 
reference.  Gas  progressing  from  left  to  right  along  the  upper  curve 
is  undergoing  cooling,  while  gas  moving  from  right  to  left  along  the 
bottom  is  undergoing  heating.  This  corresponds  to  hot  gas  cooling  as 
it  moves  away  from  the  secondary  and  cold  gas  heating  as  it 
approaches  the  secondary,  respectively. 


note,  the  small  line  dropping  out  of  the  lower  portion  of  the  curve 
is  from  the  boundary  region  of  the  problem  and  is  not  physically 
real . 

We  do  not  believe  that  this  instability  is  an  artifact  of  our 
model  alone.  Prior  to  producing  Model  1,  we  used  the  hydrodynamic 
portions  of  our  code  to  repeat  Hunt's  (1971)  problems  for  the  Mach 
1.4  and  Mach  2.4  cases.  Our  results  were  in  very  good  agreement  with 
his,  and  the  code  did  not  show  any  of  the  instability  behavior  found 
in  Model  1. 

We  have  followed  the  values  of  material  density,  energy  density 
and  temperature  in  time  within  the  instability  region  and  find  them 
to  be  very  dynamic.  Figure  5.5  displays  these  values  at  two 
different  instances.  Note  that  they  are  plotted  as  a  function  of 
zone  rather  than  distance  to  help  determine  the  source  of  the 
instability.  We  see  that  the  density  varies  inversely  with  the 
temperature,  and  that  the  energy  density  remains  relatively  smooth. 
The  onset  of  the  instability  occurs  at  about  zone  96  and  dies  out  at 
about  zone  70.  The  variations  propagate  in  a  wave  like  motion  from 
right  to  left.  Since  the  secondary  is  located  at  zone  60,  these 
variations  have  long  since  died  out  and  fluctuations  in  accretion 
rate  or  luminosity  never  appear.  We  also  see  that  the  variations  are 
spread  out  over  several  zones,  so  that  the  instability  is  not 
inherently  numerical. 

In  1977,  Hatchett  and  McCray  preformed  detailed  calculations  on 
the  transfer  of  X-rays  through  a  stellar  wind  by  assuming  either  a 
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(b).  Time  =  66730  sec. 


Figure  5.5.  Material  Density,  Temperature,  and  Energy  Density  Values 
in  the  Instability  Region.  The  density,  temperature,  and  energy 
density  are  plotted  as  a  function  of  zone  along  the  line  of  centers 
through  the  region  of  instability.  Plotted  at  two  different  times, 
the  wave  like  progression  from  right  to  left  can  be  seen.  The  minima 
and  maxima  propagate  with  the  flow  velocity. 
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constant  wind  velocity  or  a  wind  obeying  Eq.  5.6,  and  determined 
contours  of  constant  ionization  plotted  in  terms  of  the  parameter  £. 
They  found  these  contours  to  be  roughly  spherical,  with  their  centers 
displaced  from  the  source  of  the  X-rays.  For  comparison,  we  have 
produced  £,  contours  as  shown  in  Fig.  5.6.  The  shock  severely 
distorts  the  trailing  sections  of  these  contours,  showing  the  effects 
of  including  hydrodynamics.  The  right  hand  sections  are  very  nearly 
spherical,  with  the  inner  contours  more  closely  centered  about  the 
secondary  than  the  outer  ones.  The  scale  of  our  model  does  not  allow 
direct  comparison  with  their  results  beyond  the  general  agreement  in 
the  shape  of  the  contours. 

However,  we  can  use  Fig.  5.6  to  measure  the  impact  of  the  X-ray 
heating  on  the  flow.  The  £  contours  were  selected  such  that  the 
steady-state  temperatures  they  give  correspond  to  the  values  on  the 
temperature  contours  of  Fig.  5.3.  By  comparing  the  two  plots,  we 
find  that  as  we  approach  the  secondary,  the  actual  temperatures  are 
not  reached  until  much  closer  in  than  their  steady-state 
counterparts.  This  is  an  indication  that  the  X-ray  heating  effects 
are  dominated  by  the  hydrodynamic  flow. 

We  can  examine  the  impact  of  the  X-ray  heating  in  another  way. 
Let  us  compare  the  gravitational  potential  energy  with  the  thermal 
energy  of  the  gas.  If  the  thermal  energy  rapidly  approaches  the 
gravitational  energy,  capture  of  the  material  becomes  less  possible. 
To  compare  the  two  energies,  we  require  that 

GMxp/r  =  pCyT  .  (5.14) 
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Figure  5.6.  £  Contours  for  Model  1.  The  circular  shapes  of  these 
contours  outside  of  the  wake  region  are  similar  to  those  found  by 
Hatchett  and  McCray  (1977).  The  location  of  the  shock  is  best  seen 
here  where  it  lies  along  the  dips  in  the  contours  just  to  the  left  of 
the  secondary.  The  contours  along  the  upper  and  lower  edges  of  the 
plot  are  again  the  result  of  boundary  effects.  The  contour  values  in 
units  of  10“ 2 7  sec3  cm"3  are:  (1)  5.7;  (2)  3.6;  (3)  2.5;  (4)  1.7;  (5) 
1.1;  (6)  0.67. 


(5.15) 


Thus,  RQ  for  a  given  temperature  is  given  by 

Rg  -  GMx/CvT  . 

At  a  temperature  of  5  keV,  Rg  =  2.41x10^  cm,  while  at  1  keV,  Rg  is 
five  times  greater.  Re-examining  Fig.  5.3  shows  that  for  the 
spherical  regions  to  the  right  of  the  secondary,  the  1.1  keV  contour 
lies  inside  the  1  keV  Rg,  while  the  2.5  keV  contour  lies  inside  the  5 
kev  Rg.  The  drawn  out  tails  in  the  wake  pass  outside  these  radii. 
This  behavior  indicates  that  the  flow  is  dominated  by  the 
hydrodynamics  rather  than  the  X-ray  heating  in  the  region  between  the 
primary  and  the  secondary.  However,  the  X-ray  heating  begins  to 
dominate  in  the  trailing  regions  of  the  wake. 

Let  us  consider  the  time  it  takes  material  to  flow  from  Rg  of  1 
keV  to  5  keV,  and  compare  this  to  the  time  it  would  take  to  heat  the 
gas  from  1  keV  to  5  keV  over  the  same  distances.  We  developed  Eqs. 
3.19  and  3.20  to  give  us  the  X-ray  heating  time.  At  1  keV,  we  are  in 
the  Compton  dominated  heating  region,  so  we  can  use  Eq.  3.19  to 
estimate  the  heating  time  directly  from  the  Compton  time  as: 

4  Tt  r  C  T  m 

*»  •  i,  .  •  <5-16> 

Let  us  approximate  the  heating  time  by  taking  the  average  of  the 
heating  times  found  between  Rg  of  1  keV  and  5  keV.  Since  at  5  keV, 
Rg  is  2.41x10^  cm,  we  find  from  Eq.  5.16  that  the  average  heating 
time  is  1790  seconds. 

The  flow  velocity  along  the  line  of  centers  is  about  950  km 
see'*.  The  distance  between  Rg  of  1  keV  and  5  keV  is  9.64x10^  cm. 
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Thus  the  flow  crosses  this  distance  in  1014  seconds.  Along  the  line 
of  centers  at  least,  the  hydrodynamics  time  scale  is  shorter  than  the 
heating  time  by  a  factor  of  1.8.  If  the  distance  the  flow  must  cover 
increases,  by  say,  curving  in  around  behind  the  secondary,  the  flow 
time  will  increase  and  quickly  become  comparable  with  the  heating 
time.  When  this  occurs,  the  gas  will  be  able  to  heat  at  the  same 
rate  it  is  losing  gravitational  potential  energy  and  capture  will 
become  exceedingly  difficult.  From  the  temperature  structure  in  the 
wake,  we  see  that  this  is  in  fact  what  is  happening. 

The  line  accretion  model  pictures  the  bulk  of  the  material  being 
accreted  along  the  line  of  centers  from  downwind  of  the  secondary. 
The  combination  of  the  wind  force  and  the  X-ray  heating  are  acting  to 
reduce  accretion  from  this  region.  To  better  appreciate  the  flow 
near  the  secondary,  plots  of  the  velocity  field  and  mass-flux  field 
are  shown  in  Figs.  5.7  and  5.8.  Both  figures  show  a  pronounced 
reduction  in  the  flow  to  the  left  of  the  secondary,  which  corresponds 
to  the  downwind  wake  region.  Figure  5.8  was  created  by  plotting  the 
momentum  vectors  multiplied  by  their  distance  from  the  axis.  Since 
we  are  using  cylindrical  coordinates,  this  gives  us  the  contribution 
to  the  flux  integral  at  that  particular  distance.  Vectors  of 
constant  length  signify  constant  mass  flux  contributions  to  the 
accretion  rate.  In  effect,  we  are  looking  at  the  amount  of  material 
entering  each  cylindrical  shell  about  the  secondary.  Equal  vectors 
imply  that  the  same  amount  of  material  has  entered  a  given  shell. 
Thus  low  flux  over  a  large  area  contributes  as  much  material  as  a 
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Figure  5.7.  Velocity  Field  for  Model  1.  This  is  a  plot  of  the 
material  velocity  near  the  central  accreting  region  of  the  problem. 
The  lengths  are  proportional  to  a  maximum  velocity  of  2.894x10®  cm 
sec"1.  The  tail  of  each  vector  marks  the  point  at  which  the  velocity 
was  taken.  Note  the  formation  of  a  stagnation  region  near  the  rear 
of  the  secondary. 
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Figure  5.8.  Weighted  Momentum  Vectors  for  Model  1.  This  is  a  plot 
of  weighted  momentum  vectors  near  the  central  accretion  region  of 
Model  1.  The  momentum  density  vectors  are  weighted  by  their  distance 
off  axis  to  highlight  the  accretion  mass  flow.  Vectors  of  constant 
length  imply  constant  mass  flux.  Note  that  the  maximum  amount  of 
material  captured  is  from  the  sides,  rather  than  from  the  left.  All 
vectors  are  scaled  to  a  maximum  vector  of  8.68x10s  g  cm"1  sec"1. 
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large  flux  over  a  small  area.  We  note  that  most  of  the  material  is 
captured  from  the  side,  in  contradiction  with  the  line  accretion 
model,  but  in  agreement  with  the  suspected  behavior  discussed  above. 
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VI.  MODEL  TWO 


In  Chapter  5,  we  presented  a  model  for  which  the  stellar  wind 
force  and  the  X-ray  heating  combined  to  have  a  major  effect  on  the 
accretion  flow,  but  a  relatively  minor  impact  on  the  total  accretion 
rate  and  resultant  luminosity.  We  found  that  the  X-ray  heating 
effects  were  not  strong  enough  to  significantly  modify  the 
hydrodynamics  of  the  flow.  We  showed  that  the  temperatures  reached 
by  the  gas  were,  for  the  most  part,  well  below  the  corresponding 
steady  state  temperatures  of  the  gas.  We  also  found  that  the 
internal  energy  of  the  gas  was  below  the  gravitational  energy  of  the 
gas. 

A.  X-ray  Pre-heating  Analysis 

Alme  and  Wilson  (1975)  had  previously  found  that  material  flow 
could  be  significantly  affected  if  the  X-ray  heating  was  able  to 
dominate  the  hydrodynamics  of  the  flow.  The  same  results  were 
reached  analytically  by  Carlberg  (1978),  and  in  one  dimensional 
models  of  spherical  accretion  by  Cowie,  Ostriker  and  Stark  (1978). 
In  particular,  if  the  X-ray  heating  raised  the  energy  of  the  gas 
above  its  gravitational  potential  energy,  the  accretion  flow  was 
disrupted  and  the  material  was  blown  away  from  the  secondary.  We 
therefore  decided  to  create  a  second  model  which  would  show  the 
effects  of  X-ray  heating  much  more  strongly.  The  model  would  also 
correct  the  boundary  problems  which  appeared  in  Model  1. 

To  get  an  idea  of  the  parameters  we  need,  we  first  reconsider 
the  unperturbed  initial  state  of  Model  1  in  terms  of  the  steady  state 


temperatures  and  gravitational  energies.  The  radius,  Rg,  at  which 
the  gravitational  potential  energy  equals  the  thermal  energy  of  a  gas 
at  some  temperature,  T,  is  given  by  Eq.  5.18  as 

Rq  =  GMx/CvT  .  (6.1) 
Equation  4.9  can  be  rearranged  to  give  the  distance  from  the 
secondary,  Ry,  for  which  a  %  corresponding  to  a  desired  steady  state 
temperature  is  reached: 

RT  =  (£LX/  P)h  •  (6.2) 
Using  these  two  equations,  we  then  plot  the  locations  of  the  steady 
state  temperatures  and  their  corresponding  gravitational  potentials 
for  temperatures  of  1  keV  and  5  keV.  The  results  are  shown  in  Fig. 
6.1.  We  see  that  at  5  keV,  Ry  lies  just  outside  the  corresponding 
Rg,  with  Ry  averaging  1.32  Rg.  For  the  1  keV  case,  we  have  Ry 
occuring  much  farther  out  then  Rg,  averaging  2.17  Rg. 

Following  the  analysis  presented  in  Chapter  5,  we  calculate  the 
flow  time  and  X-ray  heating  time  for  the  initial  conditions.  We  find 
that  the  flow  time  is  1125  seconds,  while  the  X-ray  heating  time  is 
932  seconds.  Since  the  final  conditions  of  Model  1  have  flow  times 
less  than  the  X-ray  heating  time,  we  desire  to  construct  Model  2  such 
that  the  heating  times  dominate  the  flow  times  even  after  the  problem 
has  relaxed  to  a  steady-state.  There  appears  to  be  no  forthright  way 
of  predicting  the  values,  so  we  must  create  a  rather  extreme  set  of 
initial  conditions. 

To  determine  which  way  we  need  to  vary  the  system  parameters  to 
reach  the  desired  results,  let  us  begin  by  considering  again  Eq.  3.19 
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Figure  6.1.  Initial  Gravitational  Potential  and  Steady  State 
Temperature  Contours  for  Model  1.  Plotted  here  are  the  steady  state 
temperature  and  gravitational  potential  contours  for  the  unperturbed 
initialized  state  of  Model  1.  The  contours  were  plotted  at  1  and  5 
keV.  At  5  keV,  RT  is  1.32  times  R-,  while  at  1  keV,  RT  is  2.17  times 


ng  the  X-ray  heating  time  when  in  the  Compton  dominated  regime: 


4  it  r  m  C  T 

_ e  v 

L  <  n 
x  e 


(6.3) 


most  rapid  way  of  decreasing  the  heating  time  is  clearly  to 
ease  the  distance  to  the  secondary.  Since  these  distances  are 
g  calculated  for  various  values  of  the  gravitational  potential, 
easing  the  mass  of  the  secondary  will  have  the  desired  effect, 
naintain  self-consistency,  the  wind  velocity  near  the  secondary 
have  to  decrease,  the  wind  density  will  have  to  increase,  or  a 
ination  of  both.  At  this  point,  the  method  developed  in  Chapter 
brought  into  full  use. 


Initializing  Model  2 

We  can  use  the  results  of  Model  1  to  help  establish  the 

neters.  First,  since  the  accretion  ratio  was  0.523,  we  use  0.392 

the  accretion  efficiency  rather  than  0.75.  We  also  use  the  X-ray 

- 1 

nosity  from  Model  1,  5.212x10  ergs  sec  ,  as  the  base 

nosity.  Recalling  that  the  required  accretion  rate  for  a  given 

nosity  is  given  by  Eq.  5.1  as 

dMx/dt  =  10  Lx/c2,  (6.4) 

find  that  the  desired  luminosity  requires  an  accretion  rate  of 
16  -1 

1x10  g  sec  .  We  now  need  to  compute  the  primary  mass  loss 


with  the  aid  of  Eq.  4.11: 


dMQ/dt  = 


dMx/dt  vl  (a-Ro)2 
a  (GMx)2 


(6.5) 


start  with  the  parameters  used  in  Model  1,  and  vary  them  until  we 
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have  a  self-consistent  system.  The  values  as  finally  selected  differ 
significantly  from  the  observed  values  of  the  Vela  X-l  system.  They 
are  summarized  in  Table  6.1. 

Comparing  dMQ/dt  with  the  limit  set  by  Eq.  5.5,  we  find  that 

this  mass  loss  rate  is  a  factor  of  4  below  the  maximum  available  to 

39 

the  primary,  based  on  its  luminosity  from  Table  1  of  1.5x10  ergs 
sec-^.  We  also  note  that  the  required  primary  mass  loss  rate  lies  in 
the  middle  of  the  range  of  probable  values  established  by  the 
observations  of  Hutchings  (1976).  In  this  case  the  match  only 
demonstrates  the  reasonableness  of  the  value  since  all  of  the  other 
parameters  differ  from  the  observed.  What  we  have  established  is  an 
arbitrary,  but  self-consistent,  set  of  parameters. 

We  now  examine  the  characteristic  time  scales  given  by  these 
parameters  for  a  model  initialized  as  in  Chapter  5.  That  is,  the 
wind  is  given  a  velocity  structure  according  to  Eq.  3.5  of 

v  «  v  (i  -  R Jr)H  .  (6.6) 

oo  O 

The  secondary's  gravity  initially  has  no  effect,  and  the  gas 
temperatures  are  again  set  at  the  appropriate  steady  state 
temperatures.  The  wind  velocity  near  the  orbit  of  the  secondary  is 
427  km  sec~^,  and  Rq  at  5  keV  (re.  Eq.  6.1)  is  just  one  half  that 
found  for  Model  1,  or  1.2x10^  cm. 

With  the  above  values,  we  find  that  the  time  for  the  flow  to 
cross  from  Rg  at  1  keV  to  Rg  at  5  keV  is  1124  seconds.  The  average 
X-ray  heating  time  computed  from  Eq.  6.3  is  44.35  seconds.  The  X-ray 
heating  time  is  25  times  greater  than  the  flow  time,  so  in  fact  this 
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TABLE  6.1.  Parameters  for  Model  2 


Base  System:  4U0900-40  (Vela  X-l) 


Observed 

Model 

a  (Rg) 

52 

60 

w 

33 

41.05 

v^  (km  sec-1) 

1430 

750 

Mx  <M0> 

2 

1 

Lx  (ergs  sec-1) 

lxlO37 

5.12xl037 

Lq  (ergs  sec-1) 

1.5xl039 

3.8xl038 

dMo/dt  (Mg  yr-1) 

3.5-14.0xl0-6 

7.418xl0-1 

a 

... 

0.4 

Table  6.1.  The  parameters  used  for  Model  2  are  summarized  here.  The 
observed  values  are  taken  from  Table  2.1.  The  values  selected  for 
the  model  differ  greatly  from  the  observed  values.  The  reasons  for 
this  are  discussed  in  the  text. 


model  should  be  dominated  by  the  X-ray  heating  rather  than  the 
hydrodynamics. 

As  one  final  preliminary  analysis,  we  calculate  the  point  at 
which  the  wind  force  dominates  the  secondary's  gravity,  along  the 
line  of  centers,  and  in  the  region  of  the  wake  (re.  Eq.  5.14).  We 
find  this  point  to  be  6.13x10^  cm  from  the  secondary,  which  is  about 

4.1  R  .  This  distance  is  greater  than  that  found  in  Model  1,  but  is 

a 

closer  in  terms  of  Rg.  Considering  the  increased  scale  of  Model  2, 
the  wind  force  will  be  active  over  a  much  larger  region  of  the  wake. 

Figure  6.2  shows  a  plot  of  the  steady  state  temperatures  and  the 
gravitational  potentials  at  the  1  and  5  keV  level.  We  find  that  at  5 
keV,  the  steady  state  contour  is  pushed  to  3  Rg,  while  at  1  keV,  the 
contour  is  out  to  4.6  Rg.  Thus,  the  initial  state  configuration  is 
significantly  different  from  that  seen  in  Model  1. 

On  examining  the  velocity  structure  of  the  intial  state  we  find 
it  to  be  supersonic  inside  the  accretion  radius,  where  the  accretion 
radius  for  this  model  is  1.488x10  cm  or  nearly  a  factor  of  two 
larger  than  Model  1.  The  flow  reaches  Mach  2  at  1.03  R  ,  and  is  Mach 

a 

28.1  at  the  outer  side  boundary  which  is  53  R  .  The  velocity 

at 

structure  is  thus  similar  to  that  seen  in  Model  1. 

C.  Calculation  Results 

Model  2  is  turned  on  and  followed  in  time  until  it  too  reaches  a 
quasi -steady  state.  The  basic  results  confirm  our  suspicions  that 
the  X-rays  have  a  much  larger  effect  than  they  did  in  Model  1.  For 
the  luminosity  and  the  accretion  ratio  we  have 
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Gravitational  Potential 


—  Steady  State  Temperature 


Figure  6.2.  Initial  Gravitational  Potential  and  Steady  State 
Temperature  Contours  for  Model  2.  Shown  here  are  the  steady  state 
temperature  and  gravitational  potential  contours  for  the  unperturbed 
initialized  state  of  Model  2.  For  this  model  at  5  keV,  fL.  is  3  times 
Rq  ,  while  at  1  keV,  Ry  is  4.6  times  Rq. 
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Lx  =  1.28xl036  erg  sec-1  , 

6  =  0.235  . 

We  note  that  these  values  are  significantly  lower  than  expected,  this 
time  dropping  by  a  factor  of  5. 

With  the  interest  of  this  model  placed  on  the  X-ray  heating 
effects,  let  us  first  examine  the  temperature  contours  shown  in  Fig. 
6.3,  and  the  e,  contours  of  Fig.  6.4.  We  have  again  selected  the 
values  of  the  E,  contours  such  that  the  steady  state  temperatures  they 
corresponded  to  matched  the  temperature  contour  values.  We  first 
notice  that  the  fall  off  from  the  steady  state  levels  to  the  actual 
levels  are  not  nearly  as  drastic  as  seen  in  Model  1.  When  we  compare 
the  gravitational  potentials  with  the  temperature  contours,  we  find 
that  they  essentially  match  up  at  the  2.048  keV  level.  X-ray  heating 
is  apparently  still  much  faster  than  the  hydrodynamic  flow  times.  We 
also  see  that  the  boundary  flow  is  better  behaved,  and  that  there  is 
now  an  extended  wake  behind  the  secondary. 

We  recalculate  the  flow  and  heating  times  for  Model  2  in  its 
final  state.  The  flow  velocity  between  Rg  of  1  keV  to  5  keV  averages 
700  km  sec-*,  giving  a  flow  time  of  686  seconds.  The  corresponding 
X-ray  heating  time  calculated  from  Eq.  6.3  is  found  to  be  181.72 
seconds.  Thus  while  the  magnitude  of  the  difference  has  dropped  from 
that  of  the  initial  state,  the  X-ray  heating  time  is  still  greater 
than  the  flow  time  by  a  factor  of  3.8. 

One  effect  of  this  high  heating  rate  can  be  seen  immediately  in 
the  temperature  contours  of  Fig.  6.3.  Extremely  high  temperatures 


Figure  6.4.  £  Contours  for  Model  2.  The  C  contour  values  were  again 
selected  to  give  steady  state  temperatures  which  correspond  to  the 
values  of  the  temperature  contours  of  Fig.  6.3.  We  find  that  the 
contours  lie  much  closer  to  the  temperature  contours,  and  match  up  at 
the  2.048  keV  level  of  contour  7.  The  contour  values  in  10“26  sec3 
cm"3  are:  (1)  350.0;  (2)  7.90;  (3)  0.9;  (4)  0.45;  (5)  0.22;  (6)  0.11; 
(7)  0.043. 
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are  found  in  the  wake,  apparently  the  result  of  heated,  uncaptured 
material  flowing  past  the  secondary  and  into  the  wake  region. 
Examination  of  the  velocity  structure  of  the  gas  shows  an  interesting 
effect.  The  material  in  the  high  temperature  core  of  the  wake  is 
moving  supersonically  away  from  the  secondary  at  a  distance  of 
1.6x10**  cm.  This  is  approximately  at  the  sharp  tip  of  contour  7  in 
Fig.  6.3.  The  gas  continues  to  heat  in  the  wake,  with  the  heating 
primarily  due  to  compressional  and  viscous  effects.  The  core  of  the 
wake  reaches  velocities  of  1200  km  sec"*  and  is  supersonic  with  a 
Mach  number  of  about  2.  The  low  density  in  the  wake  keeps  it  in  the 
Compton  dominated  domain  so  that  X-ray  cooling  is  not  effective  at 
lowering  its  temperature.  The  tip  of  the  wake  apparently  is  a  second 
shock  region,  where  the  flow  changes  from  a  high  velocity,  low 
density,  high  temperature  region,  to  the  lower  velocity,  high 
density,  low  temperature  regions  of  the  surrounding  material. 

The  low  density  at  the  inner  core  of  the  wake  is  easily  seen  in 
the  density  contour  plots  of  Fig.  6.5.  Figure  6.6  again  plots  the 
energy  density  contours,  which  are  equivalent  to  pressure  contours. 
Again,  regions  of  instability  appear  as  in  Model  1.  The  cause  of  the 
onset  of  the  thermal  instability  is  more  questionable  in  Model  2  than 
it  was  in  Model  1.  The  requirement  to  extend  the  problem  mesh  leads 
to  the  creation  of  cells  with  very  high  aspect  ratios  both  on  the 
axis  and  perpendicular  to  the  secondary  at  distances  far  from  the 
secondary.  Such  cells  may  have  a  ratio  of  side  dimensions  on  the 
order  of  100.  This  causes  a  more  severe  loss  of  resolution  for  one 
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Figure  6.5.  Material  Density  Contours  for  Model  2.  Regions  of 
thermal  instability  are  again  noticable  to  the  right,  top  and  bottom 
of  the  secondary.  Their  origin  are  discussed  in  the  text.  Note  the 
formation  of  an  extended  wake  with  a  rarified  core  and  high  densities 
on  it's  boundaries.  The  contour  values  in  10”13  g  cm-3  are:  (1) 
0.25;  (2)  0.50;  (3)  1.00;  (4)  2.0;  (5)  4.0;  (6)  8.0;  (7)  16.0. 


dimension  in  those  cells  as  can  be  seen  by  the  disturbances 
introduced  into  the  energy  contours.  This  may  have  provided  the 
perturbations  needed  in  the  flow  to  start  up  the  instabilities. 

In  Fig.  6.7,  we  plot  P/F  versus  P/F  for  Model  2  as  we  did  in 
Model  1.  Again  the  essential  correctness  of  the  thermal  instability 
assumption  is  verified.  We  note  that  there  are  two  populations  of 
gas,  one  the  cool  gas  approaching  the  secondary  along  the  lower 
portion  of  the  curve,  and  the  second  the  hot  gas  which  is  moving  away 
from  the  secondary  along  the  upper  portion  of  it.  The  conspicuous 
wisp  rising  off  the  curve  along  the  107  °K  line  corresponds  to  the 
hot  thin  material  in  the  wake  behind  the  secondary.  As  noted  above, 
this  material  is  hot  gas  which  the  secondary  could  not  capture,  and 
is  expanding  away  from  the  secondary  while  undergoing  hydrodynamic 
heating.  This  is  supplying  sufficient  heating  to  keep  the  gas  hot 
and  thin  as  it  moves  away.  Remaining  thin.  X-ray  cooling  is  not 
efficient,  further  supporting  its  high  temperature.  The  high 
temperatures  and  pressures  reached  relative  to  the  X-ray  flux  place 
the  gas  in  the  wisp  seen  in  the  figure. 

This  feature  is  not  present  in  the  results  of  Model  1.  Its 
absence  may  be  due  to  the  smaller  spatial  scale  which  does  not 
include  as  large  a  segment  of  the  wake.  Looking  back  at  Fig.  5.1,  we 
see  that  the  wake  does  appear  to  be  forming  a  low  density,  high 
temperature  core,  with  a  high  density  boundary.  Expanding  Model  1 
may  lead  to  the  same  effect  that  we  see  here. 

We  did  repeat  a  time  elapsed  analysis  of  the  temperature. 


82 


p/F  (sec3  cm-3) 


Figure  6.7.  P/F  vs  p/F  for  Model  2.  This  is  a  plot  of  the  gas 
pressure  divided  by  the  X-ray  flux  versus  the  gas  density  divided  by 
the  X-ray  flux.  It  demonstrates  the  multi-valuedness  of  the 
pressure-temperature  states  of  the  gas.  Lines  of  constant 
temperature  are  included  for  reference  purposes.  The  sharpness  of 
the  corner  in  the  lower  right  of  the  plot  is  a  result  of  the 
equations  used  to  fit  the  steady  state  temperature  curves.  Gas 
moving  from  left  to  right  along  the  upper  edge  is  undergoing  cooling, 
while  gas  moving  from  right  to  left  along  the  bottom  is  gas 
undergoing  heating.  The  thin  wisp  rising  above  the  curve  along  the 
107  °K  line  represents  the  hot,  low  density  material  moving  in  the 
core  of  the  wake.  This  feature  is  not  seen  in  Model  1  since  the 
problem  did  not  include  as  much  of  the  wake. 
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density,  and  energy  density  through  the  region  of  instability  similar 
to  that  done  in  Model  1.  The  same  type  of  wave-like  behavior  was 
found  present.  The  instabilities  again  propagated  with  the  velocity 
of  the  flow. 

Figures  6.8  and  6.9  show  the  velocity  and  mass-flux  fields  for 
Model  2.  On  comparing  Fig.  6.8  with  Fig.  5.7,  we  note  that  the 
stagnation  point  behind  the  secondary  is  closer  by  a  factor  of  two. 
In  Fig.  6.9,  we  see  that  the  accretion  flow  is  again  predominately 
from  the  sides  of  the  secondary  and  there  is  essentially  no  accretion 
column  formed. 
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Figure  6.8.  Velocity  Field  for  Model  2.  This  is  a  plot  of  the 
velocity  field  near  the  central  accretion  region  of  model  2.  The 
vectors  are  scaled  to  a  maximum  vector  of  1.371x10®  cm  sec"1.  The 
tail  of  each  vector  marks  the  point  at  which  the  velocity  was  taken. 
Note  again  the  formation  of  a  stagnation  point  near  the  left  side  of 
the  secondary. 
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VII.  PRELIMINARY  ANALYSIS  OF  WIND  SHUT  DOWN  EFFECTS 

Models  1  and  2  were  created  to  examine  the  effects  of  X-ray 
heating  and  the  wind  force  on  accretion  from  a  stellar  wind.  Since 
they  did  not  allow  the  X-ray  radiation  to  impact  on  the  forces 
driving  the  wind,  the  behavior  of  the  gas  in  the  wake  is  unrealistic. 
The  densities  and  X-ray  radiation  are  such  that  the  atoms  responsible 
for  the  wind  force  are  heavily  ionized  near  the  secondary.  Thus,  we 
should  not  see  the  wind  force  dominating  the  secondary's  gravity 
through  the  major  portion  of  the  wake. 

Our  basic  goal  in  this  study  is  to  create  a  model  in  which  the 
X-rays  modify  the  forces  driving  the  wind.  With  the  behavior  of 
Models  1  and  2  established,  we  can  begin  to  formulate  the  necessary 
parameters  for  a  system  under  which  the  X-rays  ionize  the  L-shell 
electrons  responsible  for  the  wind  force.  To  do  this,  we  first  need 
to  get  some  idea  of  the  impact  that  full  wind  ionization  would  have 
on  system  modeling  requirements.  This  information  will  then  provide 
restrictions  on  our  model  parameters. 

A.  Ballistic  Particle  Model 

To  handle  this  analysis  without  explicitly  solving  the 
differential  equations  for  the  system,  we  make  the  following 
assumptions  to  create  a  worst-case  situation: 

(a).  Ionization  and  wind  turn-off  occur  on  a  time  scale 
much  smaller  than  the  time  it  takes  the  wind  to  cross 
from  the  surface  of  the  primary  to  the  orbit  of  the 
secondary.  As  a  corollary,  the  wind  force  is  assumed 
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1 


to  rapidly  turn  back  on  as  it  is  shadowed  from  the 
X-rays. 

(b) .  The  X-ray  luminosity  is  sufficient  to  strip  the 

L-shell  electrons  responsible  for  the  wind  force 
everywhere  outside  of  the  X-ray  shadow  zone. 

(c) .  The  wind  is  approximated  as  ballistic  particles  free 

to  move  under  the  wind  force  and  the  force  of 
gravity.  Hydrodynamical  effects  are  neglected. 

(d) .  The  effect  of  the  secondary's  gravity  is  ignored. 

That  these  are  worst-case  assumptions  can  be  seen  by  examining 

the  effects  of  each  assumption  in  more  detail.  For  assumption  (a), 
we  first  note  that  the  wind  crossing  times  (the  time  it  takes  for  the 
wind  to  move  from  the  primary's  surface  to  the  orbit  of  the 
secondary)  are  on  the  order  of  104  seconds,  and  the  orbital  periods 
are  on  the  order  of  10  seconds.  If  the  ionization  times  are  on  the 
order  of  the  crossing  time,  then  the  wind  will  flow  with  little 
modification  from  its  unionized  state,  and  recovery  upon  entering  the 
X-ray  shadow  will  still  occur  before  re-exposure.  Rapid  ionization 
thus  has  the  most  detrimental  effect  on  the  total  wind  flow. 

There  is  an  observational  basis  for  assumption  (a).  Conti  and 
Cowley  (1975)  performed  optical  spectroscopic  observations  of 
4U1700-37.  They  found  that  absorption  lines  were  strongly  disturbed 
in  the  region  trailing  the  secondary.  This  behavior  was  attributed 
to  the  wake  produced  by  the  secondary  as  it  moved  through  the 
primary's  wind.  This  disturbance  was  phase  dependent  and  never 
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appeared  ahead  of  the  secondary.  For  this  system,  then,  we  have  an 
actual  case  of  the  wind  flow  recovering  before  the  secondary 
reappears. 

The  effects  of  relaxing  assumption  (b)  are  more  difficult  to 
visualize.  If  the  wind  is  not  ionized  everywhere  that  it  is  exposed 
to  the  X-rays,  then  normal  wind  flow  will  occur  progressively  closer 
to  the  secondary  as  the  X-ray  luminosity  is  decreased.  One  of  the 
purposes  of  Model  3  is  to  determine  the  location  of  the  wind 
shut-down. 

As  for  assumption  (c),  inclusion  of  hydrodynamical  effects  would 
add  the  supportive  effects  of  gas  pressure,  resulting  in  a  smoothing 
of  the  flow  and  reduction  of  any  velocity  changes.  Finally,  from 
assumption  (d),  if  the  effects  of  the  secondary's  gravity  are 
included,  the  material  would  be  attracted  towards  the  secondary, 
increasing  the  likelihood  of  the  material  being  accreted  by  the 
secondary. 

The  variables  needed  for  our  ballistic  analysis  are  shown  in 
Fig.  7.1.  Region  A  is  the  X-ray  shadow  zone  and  Region  B  is  the 
X-ray  illuminated  zone.  The  binary  separation  is  given  by  a.  RQ  is 
the  primary's  radius.  R  is  the  distance  to  the  outer  edge  of  the 
shadow  zone.  4>  is  the  angle  between  the  line-of-sight  tangent  point 
and  the  line-of-centers  for  the  system.  Finally,  ¥  is  the  angle 
between  the  line-of-centers  and  R. 

When  the  wind  is  shut  down,  there  are  two  possible  conditions 
under  which  material  will  still  be  present  in  the  orbit  of  the 
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Figure  7.1.  Variables  for  Ballistic  Analysis.  Assuming  the  stellar 
wind  is  shut-off  everywhere  it  is  exposed  to  the  secondary's  X-rays, 
we  have  Region  A  as  the  shadow  zone  and  Region  B  as  the  illuminated 
zone.  a  is  the  binary  separation.  is  the  angle  to  the 
1 ine-of-sight  tangent  point  at  the  primary  surface.  ¥  is  the  angle 
at  which  the  wind  enters  Region  B  at  some  radius  R  above  from  the 
primary.  R  is  the  primary  radius. 
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secondary  for  the  secondary  to  accrete.  Let  us  define  the  crossing 

time  t  as  the  time  it  takes  the  wind  to  cross  from  the  surface  of 
c 

the  primary  at  RQ  to  the  orbit  of  the  secondary  at  a,  and  the 

intercept  time  t^  as  the  time  it  takes  the  secondary  to  move  from  ¥  = 
0  to  ¥  =  .  The  first  condition  is  then  simply  that  t  =  t..  In 
other  words,  the  secondary  directly  intercepts  material  as  it  leaves 
the  primary. 

During  the  period  t  the  wind  spends  some  time  being  driven  by 

the  wind  force  through  region  A  and  then  time  free-falling  as  the 

wind  force  is  shut  down  in  Region  B.  Therefore,  the  second  condition 

is  to  allow  the  material  to  continue  free-falling  past  the 

secondary's  orbit  to  some  radius  R^  and  then  free-fall  back  to  the 

orbit.  The  material  thus  spends  additional  time  in  free-fall  before 

crossing  the  orbit  for  the  second  time.  For  this  case,  we  also 

require  that  tc  =  t.,  but  tc  may  now  be  substantially  longer. 

It  is  clear  that  material  reaching  escape  velocity  before 

shutdown  will  never  be  captured  unless  condition  1  is  met.  Also, 

material  with  insufficient  energy  to  reach  the  secondary's  orbit  will 

not  be  captured.  This  last  situation  is  improved  if  assumption  (d) 

is  relaxed,  resulting  in  a  form  of  potential  lobe  overflow. 

With  the  limits  of  escape  velocity  and  minimum  energy 

established,  we  can  constrain  ourselves  to  considering  only  material 

leaving  between  some  v.  and  .  The  radius  at  which  the  wind 

min  max 

attains  the  minimum  energy  required  to  reach  the  secondary's  orbit 
can  be  found  from 
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(7.1) 


4  v  -  SM/Rmin  =  -GM/a 


Now,  the  wind  velocity  at  Rmin  is  specified  from  Eq.  3.5  as 

v  =  v_  (1  -  R0/R. n)h  .  (7.2) 

00  o  min' 

Equations  7.1  and  7.2  are  now  used  to  yield  an  expression  for  R  .  : 


^min 


4  vt  K  +  GM 


00  0 


4  vt  +  GM/a 


(7.3) 


Consequently,  for  we  have 

"'min  =  Arccos(R0/Rm1n)  ♦  *  .  (7.4) 

Now,  to  find  ¥  ,  we  have  the  escape  velocity  given  by 

vesc  =  (2GM/R0)?S-  (7-5) 

But  again  from  the  wind  velocity  relation  of  Eq.  7.1  we  have 

vesc  -  *-  »  *  VResc)IS  •  <7-6) 

hence, 

R«c  ■  Ro  11  -  *«c/V->'2  •  <7-7> 

Therefore,  4,mav  is  given  by 

max 

’’max  *  4rccos(RQ/Resc)  *  *  .  (7.8) 

With  V  .  and  established  to  limit  our  search  range,  we  can 
min  max 

now  proceed  to  calculate  the  required  times.  ti  is  given  simply  by 

ti  *  T  T  /  360  ,  (7.9) 

where  T  is  the  orbital  period.  Under  the  requirements  of  condition 
1,  we  first  have  the  wind  moving  under  the  wind  force  from  RQ  to  some 
radius  R.  The  time  to  cross  this  distance  is  then 

h  s  Jr0  (v»  U  "  V^V1  dr  .  (7.10) 

The  test  particles  then  free-fall  from  R  to  a,  so  this  time  is  given 
by 
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t2  =  J^(2(E+GM/r))_,s  dr  ,  (7.11) 

where 

E  =  h  vi  (1  -  R0/R)  -  GM/R  .  (7.12) 

Finally,  the  total  crossing  time  is  given  by 

tj.  =  t^  +  t2  •  (7.13) 

An  iterative  search  is  made  for  all  values  of  R  allowed  between  Y  . 

mm 

and  Ymax  attempting  to  enforce  t  =  t^.  If  no  match  is  possible,  we 
then  invoke  condition  2.  Here,  we  again  have  t^  given  by  Eq.  7.10 
but  now  t2  is  the  time  to  move  from  R  to  Rf  under  free  fall: 

t2  =  J*f(2(E  +  GM/r))'h  dr  ,  (7.14) 

with  E  given  again  by  Eq.  7.12.  We  now  need  to  include  the  free-fall 
time  from  R^.  back  to  a.  Calling  this  t^,  we  have 

t3  =  J*f  (2(GM/r  -  GM/Rf)_,s  dr  .  (7.15) 

Hence,  the  total  crossing  time  is  now 

tc  =  ti  +  t2  +  t3  .  (7.16) 

Since  we  are  able  to  compute  the  wind  velocity  from  the  above, 
it  is  possible  to  estimate  the  wind's  density  in  the  vicinity  of  the 
secondary.  Assuming  the  primary's  mass  loss  rate  is  constant,  we 
have 

dMQ/dt  =  4  tt  r*  p  v  .  (7.17) 

Thus,  at  the  point  the  wind  is  shut  down,  we  will  have  a  density  of: 


Pi  = 


dM0/dt 

2  9 
4  tt  r  v 


(7.18) 


or,  after  including  Eq.  7.1, 
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(7.19) 
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two-dimensional  model,  we  have  to  limit  the  distance  of  the  shutdown 
region  from  the  secondary.  If  the  wind  is  shut  down  too  far  away,  we 
will  not  be  able  to  adequately  account  for  the  wind  flow  prior  to  the 
intercept  of  the  secondary.  This  is  because  the  two-dimensional 
model  does  not  allow  the  wind  to  diverge  as  it  moves  radially  from 
the  primary,  so  the  lateral  extent  of  the  problem  mesh  must  remain 
small  enough  to  keep  the  flow  reasonably  planar.  For  Models  1  and  2 
we  were  not  overly  concerned  with  the  wind  behavior  far  outside  the 
accretion  radius  for  there  were  no  processes  included  in  the  models 
to  affect  the  wind  flow  along  the  side  boundaries  other  than  the 
secondary's  gravity.  We  made  the  decision  to  find  system  parameters 
which  would  cause  the  wind  shut-down  point  to  fall  between  1/4  and 
3/4  of  the  distance  between  the  secondary  and  the  primary's  surface 
for  the  initialized  state  of  the  model. 

B.  X-ray  Ionization 

We  now  need  some  way  of  determining  the  ionizing  effects  of  the 
X-ray  radiation.  With  that  determined,  we  can  estimate  the  impact  of 
the  ionization  on  the  wind  force.  The  work  of  Hatchett,  Buff  and 
McCray  (1976)  (cf.  Tarter,  Tucker  and  Salpeter  1969)  is  well  suited 
for  this  purpose.  They  examined  the  ion  population  levels  in  a  gas 
exposed  to  X-rays.  They  found  that  the  ionization  populations  could 
be  specified  as  a  function  of  the  parameter  5,  given  a  source 
spectrum  of  X-rays.  Figure  7.2  is  adapted  from  their  results  for 
oxygen  under  the  assumptions  of  an  optically  thin  gas  and  a  blackbody 
spectrum  with  a  temperature  of  2.5  keV.  The  behavior  shown  for 
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Figure  7.2.  Ionization  Equilibria  for  Oxygen.  A  sample  ion 
population  distribution  is  shown  here  for  oxygen  in  an  optically  thin 
gas  exposed  to  a  blackbody  X-ray  spectrum  with  a  color  temperature  of 
2.5  keV.  This  figure  was  adapted  from  the  results  of  Hatchett,  Buff 
and  McCray  (1976). 


oxygen  is  typical  of  the  other  elements  as  well,  although  the  exact 

features  vary.  There  is  also  a  considerable  change  in  the  behavior 

of  the  ionization  states  if  the  gas  is  optically  thick.  Oxygen,  for 

-24 

instance,  goes  from  neutral  to  fully  ionized  between  £'s  of  10  to 
10"27  if  optically  thick. 

Recall  from  Chapter  III  that  carbon,  nitrogen,  and  oxygen  are 
the  main  contributes  of  lines  responsible  for  the  wind  force 
(Castor,  Abbott  and  Klein  1975).  Our  model  code  does  not  have  the 
capability  to  explicitly  determine  the  ionization  equilibria  of  the 
elements,  so  we  made  the  approximation  that  the  wind  is  turned  off  at 
the  point  at  which  50%  of  the  oxygen  L-shell  electrons  are  stripped. 
Oxygen  is  used  primarily  because  it  is  the  only  element  for  which 
data  is  given  for  optically  thin  conditions.  Using  this  as  an 
overall  indicator  of  the  wind  ionization  state  is  not  too 
unreasonable,  for  carbon  and  nitrogen  ionize  at  slightly  larger 
values  of  S.  A  larger  source  of  error  here  is  the  fact  that  our  code 
assumes  a  22  keV  exponential  X-ray  spectrum.  Such  a  spectrum  is 
richer  in  low  energy  photons  than  a  blackbody  spectrum.  The  point  is 
that  the  model  will  show  the  main  features  of  wind  shutdown  under 
X-ray  illumination,  even  though  the  exact  details  vary. 

Let  be  the  value  of  £  for  which  50%  of  the  oxygen  L-shell 

electrons  are  stripped.  The  value  of  is  estimated  from  Fig.  7.2 
as  1 .212x10“ sec2  cm"2.  With  this  value  of  and  the  desired 
location  of  the  shutdown  point  relative  to  the  secondary  established, 
we  are  ready  to  construct  our  final  model. 
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VIII.  MODEL  3 


We  are  now  in  a  position  to  establish  the  parameters  for  Model 

3.  The  self-consistency  constraints  of  the  first  two  models  are 

still  to  hold  in  Model  3.  In  Chapter  VII,  we  discussed  the  selection 

of  the  value  of  at  which  the  wind  would  be  shut  off.  This  value 

-25  3  -3 

was  choosen  to  be  1.212x10  sec  cm  .  The  wind  shut-down  point 
must  also  fall  between  1/4  and  3/4  of  the  distance  between  the 
primary's  surface  and  the  secondary. 

A.  Initializing  Model  3 

There  is  an  additional  requirement  on  the  size  of  the  binary 
system  that  must  be  considered.  The  wind  undergoes  very  stong 
acceleration  close  to  the  surface  of  the  primary.  This  is  evident 
from  the  acceleration  found  in  the  derivation  of  the  wind  potential 
(Eq.  3.7): 

3v/3t  =  H  v^  R 0  /  r2  .  (8.1) 

The  wind  velocity  will  change  rapidly.  Consequently,  to  achieve 
reasonable  resolution  of  the  velocity,  we  have  to  have  fine  spatial 
resolution  near  the  primary's  surface.  Numerical  accuracy  in  our 
model  makes  it  desirable  to  change  cell  dimensions  no  more  than  10% 
from  one  cell  to  its  neighbor.  With  a  maximum  of  120  cells  available 
in  the  z  direction  (re.  Fig.  3.3),  there  is  a  maximum  dimension  which 
can  be  simulated  accurately.  It  turns  out  that  Vela  X-l  is  not  a 
suitable  choice  because  of  this  limitation. 

A  second  reason  exists  for  not  using  the  Vela  X-l  system.  The 
5  ff  contour  occurs  extremely  close  to  the  surface  of  the  primary 
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causing  it  to  approximate  the  worst  case  conditions  explored  in 
Chapter  VII.  It  is  clearly  apparent  that  a  new  system  is  required  as 
a  basis  for  our  model. 

The  most  obvious  candidate  to  examine  first  is  the  4U1700-37 
system.  The  parameters  for  this  system  are  given  in  Table  2.1.  It 
has  a  massive,  high  luminosity  primary  combined  with  a  small 
separation  distance  and  modest  X-ray  luminosity.  We  therefore  begin 
the  process  of  enforcing  self-consistency  using  this  system  as  a 
basis  for  model  parameter  selection. 

We  anticipate  that  the  system  will  be  initially  similiar  to 

« 

Model  1  due  to  the  higher  densities  in  the  wind.  Higher  densities 
have  the  effect  of  moving  the  contours  of  steady  state  temperatures 
closer  to  the  secondary.  From  Model  1,  we  again  assume  an  accretion 
efficiency  of  0.4.  With  the  required  accretion  rate  given  by  Eq.  5.1 
as 


dMx/dt  =  10  Lx/c2  ,  (8.2) 
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the  accretion  rate  needed  to  support  the  luminosity  of  8x10  ergs 
sec"1  (re.  Table  2.1)  is  8.89xl016  g  sec-1.  Recall  that  Eq.  4.11 
giving  the  primary  mass  loss  rate  as  a  function  of  the  accretion  rate 
had  the  form 


dMQ/dt  = 


dMx/dt  y*  (a-Ro)2 
"~~"«(GMX)* 


(8.3) 


Using  the  observed  parameters  from  Table  2.1,  Eq.  8.3  yields  a  mass 
loss  rate  of  6.7951xl0"5  Mg  yr"1. 

We  can  now  find  the  location  of  the  wind  cut-off  point.  In 
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Chapter  VII  we  selected  a  value  of  1.212x10  sec  cm  for  C0ff 
Since  rx  along  the  line  of  centers  is  simply  a-rQ,  we  can  determine  r 
with  the  aid  of  Eq.  4.10  in  the  form 


_  dMp/dt  (a-r0)g 
^  4  *  ro  Lx  v«  ^-Ro/ro)h 


(8.4) 


Doing  a  fast  iterative  search  for  r  between  Rq  and  a  locates  the 
shutdown  point  at  about  22.05  Rg,  or  just  1  Rg  above  the  surface  of 
the  primary.  This  is  entirely  too  close  for  our  model.  The  spatial 
mesh  is  not  able  to  adequately  resolve  velocities  so  near  the 
primary's  surface. 

Our  next  step  is  to  modify  the  system  parameters  in  an  attempt 
to  move  the  shutdown  point  farther  above  the  surface.  To  do  this,  we 
elect  to  reduce  Rq  from  21  Rg  to  19  Rg  without  changing  the  other 
system  parameters.  With  this  value  for  Rq,  Eq.  8.3  yields  a  required 
mass  loss  rate  of  1.0151xl0"4  Mg  yr"1.  This  mass  loss  rate  is  3.3 
times  the  maximum  mass  loss  rate  established  by  Eq.  5.4,  implying 
that  the  primary  for  this  model  must  be  3.3  times  more  luminous  than 
that  of  the  observed  system.  The  required  primary  mass  loss  rate  is 
also  a  factor  of  3.4  above  the  maximum  probable  value  from  the 
observations  of  Hutchings  (1976).  Recall  that  we  are  artifically 
enforcing  self-consistency  using  the  maximum  X-ray  luminosity.  We 
again  use  Eq.  8.4  to  determine  that  the  wind  shutdown  point  occurs  at 
about  21.8  Rg,  or  2.8  Rg  above  the  surface  of  the  primary.  This 
location  is  adequately  far  above  the  surface  of  the  primary. 
Increasing  the  distance  would  require  more  severe  modifications  of 
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the  system  parameters. 

The  minimum  spatial  resolution  obtainable  at  the  surface  of  the 
□ 

primary  is  5.0x10  cm.  This  yields  a  velocity  resolution  of  40.6  km 
sec”*.  We  can  compare  this  to  the  escape  velocity  from  the  primary's 
surface  of  775  km  sec"*  and  to  the  wind  velocity  at  the  orbit  of  the 
secondary  of  1266  km  sec"*.  Our  resolution  is  about  5%  of  the  escape 
velocity,  giving  us  good  division  between  gravitationally  bound  and 
unbound  material. 

The  parameters  finally  established  for  this  model  are  given  in 
Table  8.1.  Using  them,  we  again  initialize  Model  3  as  we  have  done 
with  Models  1  and  2.  The  accretion  radius  for  this  model  is  found  to 
be  2.48x10*°  cm.  Figure  8.1  can  be  used  to  estimate  the  initial 
effects  of  X-ray  heating.  As  one  can  see,  at  5  keV  Ry  (re.  Eqs. 
6.1-2)  actually  lies  inside  R£,  while  at  1  keV  Ry  lies  at  an  average 
distance  of  1.36  RQ.  In  this  case,  Ry  and  Rfi  show  an  even  closer 
match  than  was  found  in  Model  1  (re.  Fig.  6.1).  At  5  keV,  the  value 
of  Rg  is  1.81x10*°  cm,  and,  of  course,  at  1  keV,  Rfi  is  5  times 
larger. 

As  we  did  with  Models  1  and  2,  we  examine  the  velocity 
structures  in  the  initialized  gas.  We  find  that  the  gas  is 
supersonic  inside  the  accretion  radius,  unlike  Models  1  and  2  which 
were  subsonic  inside  the  accretion  radius.  We  can  also  calculate  the 
point  at  which  the  wind  force  exceeds  the  secondary's  gravity  along 
the  line-of-centers  on  the  downwind  side  of  the  secondary  (re.  Eq. 
5.13).  This  point  turns  out  to  be  1.88x10**  cm  from  the  secondary. 
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TABLE  8.1. 


Parameters  for  Model  3 


Base  System:  4U1700-37 


Observed 

Model 

a  (R0) 

30 

30 

W 

21 

19 

v  (km  sec'*) 

2090 

2090 

Mx  <V 

1.5 

1.5 

Lx  (ergs  sec"*) 

8xl036 

8xl036 

Lq  (ergs  sec'*) 

4xl039 

1.33xl040 

dM0/dt  (M0  yr'1) 

7.3-30.0xl0'6 

1.0l2xl0“4 

a 

.... 

0.4 

Table  8.1.  The  parameters  used  for  Model  3  are  summarized  here.  The 
observed  values  are  taken  from  Table  2.1.  The  values  selected  for 
this  model  differ  from  the  observed  values  only  in  the  radius  and 
luminosity  of  the  primary. 
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Figure  8.1.  Initial  Gravitational  Potential  and  Steady  State 
Temperature  Contours  for  Model  3.  The  steady  state  temperature  and 
gravitational  potential  contours  for  the  1  and  5  keV  levels  are 
plotted  here  for  the  initialized  state  of  Model  3.  The  5  keV 
temperature  contour  lies  inside  its  gravitational  counterpart,  while 
the  1  keV  temperature  contour  lies  an  average  of  1.36  times  the 
distance  of  its  counterpart. 


or  about  7.62  R  .  This  is  nearly  twice  as  far  in  terms  of  R  's  as 

d  a 

was  found  in  Models  1  and  2. 

B.  Calculation  Results 

We  start  the  calculation  and  again  trace  the  model's  behavior  in 

time.  Only  in  this  case,  the  model  never  approachs  the  quasi-steady 

state  of  the  previous  two  models.  Instead,  the  accretion  ratio  and 

luminosity  fluctuate  on  a  time  scale  of  about  2  hours.  The  behavior 

of  the  luminosity  is  shown  in  Fig.  8.2. 

Features  related  to  the  model  startup  period  cover  the  interval 

from  0  hours  to  about  13  hours.  This  is  the  time  needed  for  the 

initial  material  to  clear  the  central  portion  of  the  spatial  mesh. 

After  the  startup  period,  we  note  that  sharp  peaks  occur  in  the 

38  -1 

luminosity.  The  maximum  luminosities  are  4  to  5x10  ergs  sec  , 

36  -1 

while  the  minimum  luminosity  is  8x10  ergs  sec  . 

Referring  back  to  Eq.  2.7,  we  find  that  Ledd  for  this  model  is 
38  1 

2x10  ergs  sec  .  Examining  the  duration  of  the  peaks,  we  note  that 
they  exceed  Lgdd  for  periods  of  about  22  minutes.  We  will  address 
this  problem  of  exceeding  the  Eddington  luminosity  in  Chapter  IX. 
However,  we  note  here  that  the  model  did  not  incorporate  a  treatment 

of  radiation  pressure,  hence  it  could  not  be  expected  to  properly 

account  for  the  limiting  effects  of  the  Eddington  luminosity. 

Because  of  the  variations  occuring  in  the  luminosity,  the 
contour  and  vector  plots  for  this  model  represent  only  one  instant  in 
time.  All  of  these  plots  are  made  at  24.29  hours  into  the  model  run. 
This  time  represents  a  period  of  decreasing  luminosity  immediately 
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following  the  last  peak  in  Fig.  8.2.  It  also  corresponds  to  the  time 
at  which  we  terminated  our  calculation.  To  examine  the  time 
variations,  we  produced  several  movies  of  the  density  contours 
evolving  in  time.  Analysis  of  these  movies  is  used  in  the 
interpretation  of  the  plots  presented  here. 

The  density  contour  plot,  shown  in  Fig.  8.3,  exhibits  very 
complicated  structure.  The  closing  of  the  contours  to  the  right  of 
the  secondary  is  a  direct  result  of  the  action  of  the  varying  X-ray 
luminosity  on  the  wind  force.  A  feedback  mechanism  is  evidently  in 
action.  The  varying  density  in  the  wind  as  it  passes  the  secondary 
causes  a  variation  in  the  luminosity.  Increasing  luminosity  causes 
the  wind  shutdown  point  to  move  closer  to  the  primary's  surface,  and 
vice  versa.  This  in  turn  introduces  new  fluctuations  in  the  wind 
density. 

To  better  visualize  the  enhancement  process,  let  us  follow  a 
packet  of  wind  material  as  it  travels  from  the  surface  of  the  primary 
to  the  secondary.  As  this  packet  leaves  the  surface  of  the  primary, 
it  reaches  a  state  at  which  5  is  less  than  At  this  point  it 
begins  to  decelerate  under  the  influence  of  the  primary's  gravity. 
Other  packets  are  able  to  catch  up  to  it  from  behind,  causing  the 
density  of  the  packet  to  increase.  As  a  consequence  5  now  increases 
above  and  the  packet  can  move  again  under  the  wind  force  until 
the  force  is  shut  off  again. 

This  effect  causes  waves  of  enhanced  density  to  be  created.  As 
these  sweep  by  the  secondary,  they  provide  additional  material  for 
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Figure  8.3.  Density  Contours  for  Model  3.  As  in  the  previous 
models,  the  secondary  is  marked  by  the  small  x.  The  surface  of  the 
primary  now  lies  essentially  on  the  right  hand  border.  The  material 
flows  from  right  to  left.  Model  3  exhibited  strong  oscillations,  and 
thus  the  contours  changed  with  time.  Those  presented  for  Model  3 
were  taken  at  24.29  hours  of  simulation  time  into  the  problem.  We 
note  in  particular  the  build  up  of  enhanced  density  regions  between 
the  secondary  and  the  right  hand  side,  and  the  strong  turbulence 
introduced  in  the  wake.  These  enhanced  density  regions  are  both 
caused  by,  and  result  from,  the  fluctuating  luminosity.  The  contour 
values  in  units  of  10'11  g  cm"3  are:  (1)  0.1;  (2)  0.22;  (3)  0.46;  (4) 
1.0;  (5)  2.2;  (6)  4.6;  (7)  10.0 
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accretion,  causing  an  increase  in  the  luminosity.  The  increased 
luminosity  in  turn  affects  the  location  where  the  wind  is  shut  off. 
This  feed-back  mechanism  causes  the  luminosity  to  fluctuate  on  the 
same  time  scale  as  the  wind-crossing  time.  The  effects  seen  in  an 
actual  system  will  be  modified  by  the  effects  of  orbital  motion  and 
any  material  storage  time  introduced  by  formation  of  an  accretion 
disk  around  the  secondary.  Since  these  enhanced  density  waves  are 
propelled  by  the  wind  potential,  it  is  not  suprising  that  the  major 
flaring  is  occurring  on  the  same  time  scale  as  the  wind  crossing 
time. 

The  energy  density  contours  are  displayed  in  Fig.  8.4.  When 
followed  in  time,  the  energy  density  contours  display  the  behavior 
expected  in  light  of  the  fluctuations  observed.  As  the  luminosity 
increases,  the  energy  density,  hence  the  pressure,  is  seen  to 
increase.  This  inhibits  the  flow  of  incoming  material,  causing  a 
drop  in  the  accretion  rate.  As  the  accretion  flow  is  shut  off,  the 
luminosity  drops,  the  gas  cools,  and  the  accretion  rate  picks  up 
again. 

Figures  8.5  and  8.6  display  the  temperature  and  5  contours.  As 
in  the  previous  two  models,  the  values  of  the  S  contours  were 
selected  such  that  the  steady  state  temperatures  they  represent 
correspond  to  the  values  of  the  temperature  contours.  For  the  time 
shown,  we  find  that  the  two  sets  of  contours  essentially  match  up, 
except  near  the  inner  contours  which  represent  the  maximum 
temperatures.  On  examining  plots  made  during  a  period  of  increasing 
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Figure  8.4.  Energy  Density  Contours  for  Model  3.  The  contour  values 
in  units  of  ergs  cm"3  are:  (1)  7.68;  (2)  38.4;  (3)  192.0;  (4)  960.0; 
(5)  4800.0;  (6)  24000.0 


Figure  8.5.  Temperature  Contours  for  Model  3.  Strong  turbulence  is 
clearly  seen  in  the  wake  behind  the  secondary.  Note  the  radical 
departure  from  the  smooth  spherical  shapes  of  the  previous  models. 
The  contour  values  in  eV  are:  (1)  1.7;  (2)  5.1;  (3)  15.0;  (4)  46.0; 
(5)  140.0;  (6)  410.0;  (7)  1200.0 


luminosity,  we  find  that  the  steady  state  temperatures  fall  well 
within  the  actual  temperatures  of  the  gas.  This  indicates  that  X-ray 
heating  dominates  the  flow  as  the  luminosity  approaches  its  peak 
values. 

Figure  8.7  is  a  P/F  versus  p/F  plot.  We  again  find  two  gas 

populations  present.  The  extent  of  the  problem  mesh  provides  a 

larger  range  of  p/F  values  than  was  found  in  the  earlier  models.  We 

note  that  the  gas  is  well  behaved  in  its  very  hot  and  cool  states, 

but  that  the  two  gas  populations  are  again  limited  to  the  region 
-27  -25 

between  p/F  =  10  and  p/F  =  10  .  Tracing  this  plot  in  time  shows 

an  interesting  behavior.  As  the  luminosity  begins  to  rise  after  a 
minimum,  the  distinction  between  the  two  populations  becomes  blurred. 
However,  as  the  luminosity  approaches  its  peak,  material  begins  to 
jump  sharply  upwards  until  the  two  gas  populations  are 
re-established.  This  creation  of  two  populations  within  the  gas  is 
again  the  cause  of  the  instability  regions  seen  in  Figs.  8.4-6  lying 
to  the  right  of  the  secondary. 

Finally,  in  Figs.  8.8  and  8.9,  we  display  the  velocity  field  and 
weighted  momentum  vectors  near  the  secondary.  Taking  into  account 
the  fact  that  the  luminosity  has  just  passed  a  peak,  we  are  not 
surprised  to  see  the  extended  stagnation  regions  behind  the 
secondary.  The  scale  of  these  plots  is  essentially  the  same  as  that 
used  for  Models  1  and  2  (re.  Figs.  5.7,  5.8,  6.8,  and  6.9).  While 
difficult  to  see  on  the  velocity  plots,  the  momentum  plot  (Fig.  8.9), 
shows  a  pronounced  outward  flow  of  material.  This  was  not  observed 
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Figure  8.7.  P/F  vs  p/F  for  Model  3.  As  with  Models  1  and  2,  the 
presence  of  two  gas  populations  is  evident  by  the  separation  of  the 
curve  near  it's  hump.  The  sharp  corner  in  the  lower  right  is  an 
artifact  of  the  fit  to  the  steady  temperature  curve  (Fig.  3.1). 


Figure  8.8.  Velocity  Field  for  Model  3.  The  vectors  are  drawn 
proportional  to  the  maximum  vector  of  2.058x10s  cm  sec-1.  The 
velocities  are  those  computed  at  the  location  of  the  tail  of  the 
vectors. 
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Figure  8.9.  Weighted  Momentum  Vectors  for  Model  3.  The  momentum 
vectors  indicate  the  momentum  at  the  location  of  the  tails  multiplied 
by  the  distance  off  axis.  This  serves  to  emphasize  the  mass  flux 
into  the  accretion  region.  Vectors  of  equal  length  signify  equal 
mass  flow.  Note  that  most  of  the  accreted  material  is  captured  from 
the  sides  in  agreement  with  the  line  accretion  model.  The  vectors 
are  scaled  to  a  maximum  vector  of  6.942xl05  g  cm-1  sec-1. 
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in  Models  1  and  2.  The  outward  flow  material  to  the  left  of  the 


secondary  could  be  responsible  for  the  clumping  of  material  off  the 
line  of  centers  to  the  left  of  the  secondary  as  see  in  Fig.  8.3. 

The  most  recent  series  of  observations  on  the  4U1700-37  system 
were  presented  by  Pietsch  (1980),  and  Dolan  (1980).  Pietsch  found 
that  the  4U1700-37  system  exhibits  flares  in  its  luminosity  of 
magnitude  on  the  order  of  25,  with  time  scales  varying  from  minutes 

OC  1 

to  hours.  The  maximum  luminosity  observed  was  5.0x10  erg  sec  . 

We  compare  this  to  the  flares  exhibited  by  Model  3  which  vary  in 

luminosity  by  factors  of  30  to  100,  and  on  similar  time  scales. 

38 

Although  the  peak  luminosity  after  the  start-up  period  is  5.0x10 
ergs  sec"^,  this  is  related  to  the  fact  that  the  maximum  luminosity 
was  used  to  establish  the  model  parameter  set.  In  the  discussion 
which  follows  in  Chapter  IX,  we  will  address  how  Model  3  can  be 
adjusted  to  bring  it  closer  in  line  with  observations. 


IX.  DISCUSSION 

The  purpose  of  this  study  was  to  combine  recent  results  in  three 
areas  pertaining  to  X-ray  binary  star  systems:  first,  hydrodynamic 
studies  of  mass  exchange  mechanisms;  second.  X-ray  radiation 
transport  through  optically  thin  gases;  and  third,  new  models  for 
radiation  driven  stellar  winds  in  early  type  supergiant  stars.  In 
Chapter  III,  we  presented  methods  for  incorporating  X-ray  heating 
effects  and  the  stellar  wind  force  into  the  hydrodynamic  equations. 

Models  1  and  2  were  presented  in  Chapters  V  and  VI.  Both  models 
included  X-ray  heating  and  the  stellar  wind  force.  However,  in 
neither  model  was  ionization  of  the  wind  by  the  X-rays  allowed  to 
modify  the  strength  of  the  wind  force.  From  these  models  we  have 
learned  that  X-ray  heating  and  the  wind  force  can  cause  reductions  in 
the  accretion  rates  and  resulting  luminosities  as  compared  to 
predictions  of  previous  numerical  hydrodynamic  models  and  Hoyle  and 
Lyttleton's  (1939)  analytic  line  accretion  model.  The  anticipated 
impact  of  significant  X-ray  pre-heating  of  the  gas  was  specifically 
verified  in  Model  2. 

In  Chapter  VII  we  developed  a  simplified  model  which 
demonstrated  some  of  the  characteristics  of  a  stellar  wind  flowing 
under  the  influence  of  intense  X-ray  radiation.  This  led  to  Model  3 
which  was  discussed  in  Chapter  VIII.  Model  3  represents  a  step 
beyond  previous  calculations  in  that  the  strength  of  the  wind  force 
was  allowed  to  be  modified  by  X-ray  ionization.  The  work  of 
Hatchett,  Buff  and  McCray  (1976)  allowed  us  to  estimate  the 


ionization  state  of  the  wind.  They  determined  the  ionization  state 
of  a  gas  of  cosmic  abundance  as  a  function  of  the  same  variable, 
which  we  used  for  our  heating  parameter.  This  enabled  us  to  estimate 
the  point  at  which  the  wind  force  was  shut  off.  The  substantial 

effects  on  the  flow  led  to  periodic  flaring  of  the  X-ray  luminosity. 

The  shapes  of  the  ionization  regions  were  also  strongly  modified  as  a 
result  of  the  varying  wind  densities.  The  peak  luminosities 

resulting  from  Model  3  are  artifically  high.  The  model  was 
formulated  to  be  self-consistent  using  the  peak  luminosity.  The 

unexpected  strong  variations  in  the  wind  density  led  to  the  large 

luminosity  fluctuations. 


A.  Model  Scaling  Analysis 

At  this  point  we  can  consider  how  the  results  of  Model  3  would 
scale  with  changes  in  system  parameters  under  steady  state 
conditions.  We  first  note  that  the  wind  density  is  directly 
proportional  to  the  mass  loss  rate  from  Eq.  4.7: 


P  » 


dM0/dt 


(9.1) 


4  tt  r  v 

Since  the  accretion  rate  is  proportional  to  the  wind  density  from  Eq. 


4.3: 


dMx/dt  =  11  Rg  px  vx  a  '  0-2) 
the  X-ray  luminosity  can  be  expected  to  be  proportional  to  the  mass 
loss  rate.  Additionally,  the  heating  parameter,  E,  is  directly 
proportional  to  the  density,  and  inversely  proportional  to  the  X-ray 
luminosity  from  Eq.  3.10: 
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K  =  P  /  Lx  .  (9.3) 
Thus,  £  would  remain  constant  at  any  given  point  as  the  mass  loss 
rate  was  varied.  Since  £  is  constant,  we  expect  no  differences  in 
the  wind  ionization  structure  as  the  mass  loss  varies.  The  steady 
state  temperatures  in  the  wind  will  also  remain  unchanged.  However, 
since  the  X-ray  luminosity  enters  into  the  calculation  of  the  X-ray 
heating  time,  we  need  to  examine  the  impact  of  a  varying  mass  loss 
rate  on  the  X-ray  heating  rate  in  more  detail. 

To  estimate  the  changes  in  the  X-ray  heating  rate,  recall  that 
the  X-ray  heating  time  was  given  by  Eq.  3.20  as 


tx  = 


1  +  (Et/E)‘ 


(9.4) 


where  the  Compton  heating  time  is  given  by  Eq.  3.19  as 


tc  = 


4  n  rc  me  Cy  T 
Lx  " 


(9.5) 


The  rate  of  change  in  energy  due  to  X-ray  heating  is  given  by  Eq. 
3.21  as 


(3E/3t)r  =  (Ess  -  E)/tx  .  (9.6) 
We  see  from  Eqs.  9.4  and  9.5  that  the  X-ray  heating  time  is  inversely 
proportional  to  the  X-ray  luminosity,  and,  therefore  also  inversely 
proportional  to  the  primary's  mass  loss  rate.  If  the  X-ray  heating 
time  increases,  Eq.  9.6  indicates  a  corresponding  decrease  in  the 
rate  of  change  in  energy.  Apparently,  decreasing  the  mass  loss  rate 
while  maintaining  self-consistency  in  the  luminosity,  will  have  the 
effect  of  decreasing  the  impact  of  X-ray  heating. 
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The  above  analysis  has  been  applied  to  steady  state  conditions. 
To  extend  this  analysis  to  dynamic  conditions  such  as  found  in  Model 
3,  we  begin  by  noting  that  if  S  remains  constant  at  any  point,  the 
same  flare  producing  process  will  be  present.  The  wind  force 
shut-down  will  occur  at  the  same  locations  in  the  flow,  leading  to 
similar  regions  of  material  density  enhancements.  This  leads  us  to 
conclude  that  flaring  of  Model  3  should  continue  to  be  present  with 
the  same  general  features  as  the  primary's  mass  loss  is  decreased. 

Being  able  to  scale  the  luminosity  with  the  primary's  mass  loss 
rate  allows  us  to  consider  the  effects  of  such  a  reduction  on  Model 
3.  If  we  reduce  the  mass  loss  rate  by  a  factor  of  100,  the  peak 
luminosities  come  into  agreement  with  the  observations  of  the  base 
system,  4U1700-37.  The  duration  of  the  flares,  and  the  behavior  in 
the  wake,  which  are  dominated  by  X-ray  heating,  can  not  be 
determined.  A  reduction  by  a  factor  of  100  in  the  primary's  mass 
loss  rate  results  in  a  rate  that  is  a  factor  of  30  below  the  maximum 
mass  loss  rate  available  to  the  primary. 

The  lowering  of  the  maximum  luminosity  also  alleviates  the 
problem  of  exceeding  the  Eddington  luminosity,  Lgdd.  After  the 
initial  start  up  period.  Model  3  exceeds  Ledd  by  factors  of  about  3 
on  three  occasions  for  periods  of  22  minutes.  Since  radiation 
pressure  is  not  accounted  for  in  our  hydrodynamic  equations.  Model  3 
can  not  properly  limit  accretion  as  Ledd  is  reached.  However,  the 
durations  of  the  X-ray  luminosity  peaks  are  short  enough  that  they 
may  be  observed  even  if  radiation  pressure  is  included.  This  will 
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depend  on  how  smooth  the  flow  is  by  the  time  it  reaches  thfr  neutron 
star's  surface.  If  the  flow  is  greatly  disrupted  from  spherically 
symmetric  flow,  brief  outbursts  exceeding  L  ^  may  occur  as  was 
discussed  by  Vitello  (1978). 

B.  Observational  Comparisons  with  Model  3 

The  most  recent  observational  work  on  4U1700-37  was  done  by 
Pietsch  et  al.  (1980)  and  is  summarized  in  Table  2.1.  They  found 
slowly  varying  X-ray  flaring  time  scales  of  about  0.5  to  1.0  hour. 
This  flaring  did  not  appear  to  be  periodic  in  time.  They  also  found 
that  the  X-ray  spectrum  could  be  fit  by  a  28  keV  thermal 
bremsstrahlung  spectrum  in  the  20-180  keV  energy  range. 

Model  3  predicted  periodic  X-ray  flaring  on  a  time  scale  of 
about  3.5  hours.  Since  this  period  is  about  equal  to  the  wind 
crossing  time,  it  may  be  related  to  the  way  the  wind  force  is 
disrupted.  The  high  luminosity  peaks  are  able  to  turn  the  wind  force 
off  very  close  to  the  primary's  surface.  This  requires  the  regions 
of  enhanced  density  to  travel  the  full  distance  from  the  primary  to 
the  secondary.  A  more  accurate  model  of  the  wind  ionization 
structure  would  allow  the  wind  force  to  turn  off  gradually.  The 
density  enhancements  would  then  form  closer  to  the  secondary  and  have 
higher  velocities  than  our  model  predicts.  The  regions  would  also 
tend  to  have  longer  density  scale  lengths,  making  the  actual  impact 
on  the  flaring  behavior  difficult  to  predict. 

The  observed  spectrum  is  harder  than  the  22  keV  exponential 
spectrum  which  was  assumed  for  our  steady-state  temperature  curve 
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(Fig.  3.1).  This  would  cause  the  wind  to  be  ionized  farther  from  the 
secondary  than  in  our  model.  Again  the  character  of  the  density 
enhancements  would  be  changed. 

C.  Weaknesses  of  Model  3 

There  are  three  major  weaknesses  in  Model  3.  The  first  is  the 
restriction  of  the  model  to  two  dimensions,  as  was  discussed  in 
Chapter  III.  The  second  is  the  lack  of  any  orbital  motion.  The 
symmetry  invoked  for  the  two  dimensional  models  made  it  impractical 
to  include  such  motion.  This  latter  omission  impacts  most  heavily  on 
the  dynamics  of  the  wake,  while  the  features  found  near  the  accretion 
region  are  not  likely  to  be  seriously  affected.  The  ballistic  model 
of  Chapter  VII  addressed  the  three  dimensional  nature  of  the  actual 
wind  flow.  It  is  apparent  that  a  three  dimensional  model  capable  of 
handling  the  extended  wind  structure,  along  with  the  orbital  motion, 
would  be  required  for  a  complete  modeling  of  the  processes  in  X-ray 
binary  systems. 

The  final  weakness  lies  in  our  method  of  selecting  a 
representative  value  of  our  heating  parameter,  at  which  the 
wind  force  is  assumed  to  be  shut  off.  As  noted  in  Chapter  VIII,  £ 
covers  a  large  range  of  ionizations  states,  going  from  neutral  to 
fully  ionized  as  £  varies  by  about  1000.  Hence,  the  wind  force  is 
gradually  weakened  rather  than  being  switched  off  suddenly.  In 
particular,  estimates  made  from  the  results  of  Hatchett  and  McCray 
(1977)  indicate  that  L-shell  ionization  would  begin  at  the  surface  of 
the  primary  and  be  complete  near  the  secondary. 
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The  effect  of  suddenly  switching  off  the  wind  force  is  to  create 
a  sharp  density  variation  in  the  wind  which  passes  the  secondary  at  a 
velocity  lower  than  the  normal  wind  velocity.  Gradually  weakening 
the  wind  force  would  produce  a  smoother  density  variation  in  the 
wind,  but  it  would  be  moving  at  a  higher  velocity.  These  two 
processes  can  produce  the  same  flaring  characteristics.  The  flaring 
behavior  is  dependent  on  the  maximum  density  reached  in  the  enhanced 
density  region  and  on  the  time  it  takes  this  region  to  pass  the 

secondary.  In  other  words,  the  flaring  is  the  same  if  the  density 

pulse  seen  by  the  secondary  is  the  same. 

D.  Validity  of  the  Optically  Thin  Assumption 

It  is  now  appropriate  to  discuss  the  validity  of  our  assumption 
that  the  wind  flows  are  optically  thin.  We  can  check  this  assumption 
by  calculating  the  optical  depth  for  Compton  scattering  along  3  paths 
through  the  models.  We  can  then  compute  an  adjustment  factor  to  give 
the  total  opacity  as  a  function  of  frequency  relative  to  the  Compton 
optical  depth.  The  Compton  optical  depth  gives  an  indication  as  to 

the  validity  of  the  optically  thin  assumption  and  the  accuracy  of  the 

heating  model  in  various  regions  of  the  problem  mesh.  The  adjusted 
depth  determines  the  degree  to  which  the  X-ray  spectrum  is  modified 
from  the  assumed  22  keV  exponential  shape.  This  will  allow  a 
judgement  to  be  made  concerning  the  accuracy  of  the  ionization 
structure  predicted  for  the  wind. 

We  estimate  the  Compton  optical  depth  using  the  relation 

x  =  <e  /  p  dx  ,  (9.7) 
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where  Kg  is  the  Thompson  electron  scattering  opacity.  For  cosmic 

2  -1 

abundances,  <Q  is  0.33  cm  g  .  To  estimate  the  frequency  dependence 
of  the  opacity,  we  use  the  analytical  fits  given  by  Basko  and  Sunyaev 
(1973).  These  fits  smooth  over  the  absorption  edges  of  the  L  and  K 
shells,  but  are  satisfactory  for  our  purposes.  Basko  and  Sunyaev 
found  that  for  photon  energies  below  0.87  keV,  the  opacity  was  given 
by 

■c(v)  =  7.3xl023  (vQ/v)3  Nq  A"1  ,  (9.8) 

where  Nq  is  Advogadro's  number,  A  is  the  atomic  weight  of  the 
material,  and  vQ  is  equal  to  1  keV.  For  cosmic  abundances,  A  is 
equal  to  1.16198,  so  Eq.  9.7  can  be  expressed  as 

k(v)  =  114.65  (vQ/v)3  Ke  .  (9.9) 

For  photon  energies  above  0.87  keV,  Basko  and  Sunyaev  give 

k(v)  =  2.5xl0’22  (v0/v)2*4  Nq  A-1  ,  (9.10) 

or, 

k(v)  =  130  (v0/v)2,4  <e  .  (9.11) 

Equations  9.9  and  9.11  are  plotted  in  Fig.  9.1. 

The  Compton  optical  depths  are  computed  for  each  model  along 
three  different  paths.  Path  1  is  from  the  secondary  to  the  primary 
along  the  1 ine-of-centers .  Path  2  is  from  the  secondary  along  a  path 
perpendicular  to  the  1 ine-of-centers.  Finally,  path  3  is  away  from 
the  primary  along  the  1 ine-of-centers,  moving  out  through  the  wake 
* ogion.  The  resulting  optical  depths  are  plotted  in  Figs.  9.2,  9.3, 
and  9.4  for  Models  1,  2  and  3,  respectively. 

On  examining  Fig.  9.2,  we  see  that  Model  1  is  very  thin  to 
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Figure  9.1.  Optical  Depth  Factor  vs.  Photon  Energy.  This  graph 
provides  a  means  for  estimating  the  frequency  dependence  of  the 
optical  depth.  The  Optical  Depth  Factor  represents  the  increase  in 
depth  a  photon  of  given  energy  experiences  over  the  Compton 
Scattering  Depth. 
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Figure  9.2.  Compton  Optical  Depths  for  Model  1.  Compton  optical 
depths  are  plotted  along  three  paths  with  the  secondary  as  the 
starting  point.  Path  1  is  along  the  line-of-centers  towards  the 
secondary.  Path  2  is  perpendicular  to  the  line-of-centers.  Path 
three  is  along  the  line-of-centers  away  from  the  primary.  The  arrow 
marks  the  distance  at  which  the  material  temperature  along  each  path 
drops  below  0.1  keV. 
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Figure  9.3.  Compton  Optical  Depths  for  Model  2.  Compton  optical 
depths  are  plotted  along  three  paths  with  the  secondary  as  the 
starting  point.  Path  1  is  along  the  line-of-centers  towards  the 
secondary.  Path  2  is  perpendicular  to  the  line-of-centers.  Path 
three  is  along  the  line-of-centers  away  from  the  primary.  The  arrows 
mark  the  points  along  the  paths  at  which  the  material  temperatures 
drops  below  0.1  keV. 
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Figure  9.4.  Compton  Optical  Depths  for  Model  3.  Compton  optical 
depths  are  plotted  along  three  paths  with  the  secondary  as  the 
starting  point.  Path  1  is  along  the  line-of-centers  towards  the 
secondary.  Path  2  is  perpendicular  to  the  line-of-centers.  Path 
three  is  along  the  line-of-centers  away  from  the  primary.  The  arrow 
marks  the  point  along  each  path  at  which  the  material  temperature 
drops  below  0.1  keV. 


Compton  scattering.  The  maximum  depth  is  0.033  Compton 
mean-free-paths.  Referring  back  to  Fig.  9.1,  we  see  that  the 
increase  in  optical  depth  by  the  factor  of  100  needed  to  make  the 
material  optically  thick  is  obtained  for  photons  with  energies  of  1 
keV  or  less.  Thus  our  optically  thin  assumption  holds  reasonably 
well  in  all  regions  of  Model  1. 

The  Model  2  results  are  shown  in  Fig.  9.3.  We  find  that  the 

model  becomes  optically  thick  along  the  paths  1  and  2  for  distances 

greater  than  about  10^  cm.  There  is  a  drastic  change  in  the  optical 

10 

depth  of  the  wake,  as  indicated  by  the  sudden  jump  occurring  at  10 

cm  along  path  3.  This  corresponds  to  the  very  tip  of  the  hot  central 

wake.  Beyond  this  point  the  cool,  dense  material  has  finally 
collapsed  onto  the  axis,  causing  the  very  rapid  increase  in  optical 
depth  from  that  point  outward.  This  behavior  is  limited  to  the  dense 
material  on  the  axis.  Material  outside  the  wake  has  densities 
similar  to  that  found  along  the  path  2. 

The  far  wake  is  optically  thick  not  only  to  Compton  scattering, 
but  also  to  photons  with  energy  well  above  10  keV.  The  material 
along  path  3  shows  a  Compton  depth  of  0.32  mean-free-paths.  From 
Fig.  9.1,  we  see  that  this  material  becomes  optically  thick  to 
photons  with  energies  of  4.5  keV  or  less.  Since  the  densities  in 

Model  2  are  not  that  different  from  those  of  Model  1,  Model  1  would 

have  the  same  cut-off.  Model  1  therefore  predicts  a  spectral  cut-off 
very  similiar  to  the  2. 2-4. 4  keV  cut-off  actually  observed  in  the 
Vela  X-l  system  (re.  Table  2.1). 
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To  judge  the  impact  on  the  heating  model  caused  by  the  high 
optical  depths,  we  have  indicated  on  Fig.  9.3  the  distances  along 
each  path  at  which  the  temperatures  drop  below  0.1  keV.  At  the 
distances  marked,  the  steady  state  temperatures  are  very  low,  so  that 
very  little  X-ray  heating  occurs.  Thus,  X-ray  heating  does  not 
significantly  effect  the  flow  beyond  the  region  of  the  model  where 
the  optically  thin  assumption  holds. 

Finally,  Fig.  9.4  shows  the  behavior  of  the  optical  depths  in 
Model  3.  The  behavior  is  nearly  identical  to  that  seen  in  Model  2, 
except  that  the  higher  densities  encountered  in  Model  3  have  shifted 
the  features  upward  by  an  order  of  magnitude.  We  again  see  a  sharp 
jump  in  the  optical  depth  along  the  wake  (path  3),  while  the  other 
two  paths  show  a  smooth  increase.  We  have  also  marked  with  an  arrow 
the  distance  at  which  the  temperatures  along  the  three  paths  drop 
below  0.1  keV.  At  this  point,  the  flow  has  become  optically  thick  to 
photons  of  about  3  keV.  Thus,  our  optically  thin  assumption  breaks 
down  much  more  rapidly  for  Model  3  than  for  the  first  two  models. 
The  region  of  strong  X-ray  heating  near  the  secondary  is  still 
reasonably  optically  thin.  However,  beyond  a  distance  of  2x10^  cm 
the  flow  is  decidedly  optically  thick.  This  implies  that  the 
ionization  structures  given  by  Model  3  are  not  entirely  accurate  in 
this  region.  In  particular,  the  wind  shut-down  region  is  found  to 
move  into  and  out  of  the  optically  thin  region  as  the  X-ray 
luminosity  varies. 

Utilizing  Fig.  9.1,  we  find  a  spectral  cut-off  of  about  7.5  keV 


131 


for  regions  lying  outside  of  the  wake.  A  2. 2-5. 5  keV  cut-off  was 
noted  in  Table  2.1  for  the  4U1700-37  system.  This  is  another 
indication  that  the  densities  in  Model  3  are  higher  than  the 
densities  of  the  base  system.  This  is  not  suprising,  since  the 
primary's  mass  loss  rate  was  artificially  high  to  enforce 
self-consistency  at  the  maximum  observed  luminosity.  Hence  the 
artificially  high  wind  densities  are  expected. 

Some  care  must  be  exercised  in  comparing  the  optical  depths  of 
these  models  with  the  real  systems.  In  the  actual  systems,  the  wind 
flow  is  spherically  symmetric  with  respect  to  the  primary.  Our 
models  used  planar  flows  which  do  not  exhibit  radial  divergence. 
Thus,  the  wind  density  in  our  models  remains  constant  along  paths 
perpendicular  to  the  line  of  centers.  Spherical  flow  would  yield 
continually  decreasing  densities  along  similar  paths.  This  causes 
our  models  to  predict  optical  depths  higher  than  those  expected  in 
the  actual  star  systems. 

While  all  three  models  show  a  tendency  to  become  optically  thick 
far  from  the  secondary,  we  note  that  the  temperature  of  the  material, 
where  it  becomes  optically  thick,  is  on  the  order  of  a  few  eV.  Thus, 
the  regions  near  the  secondary  where  the  flow  is  most  strongly  heated 
by  the  X-rays  are  still  optically  thin,  so  that  our  X-ray  heating 
model  holds  rather  well.  However,  the  description  of  the  ionization 
structure  of  the  gas  beyond  this  region  becomes  progressively  more 
inaccurate  with  increasing  optical  depth.  This  is  due  to  absorption 
of  the  lower  energy  photons  which  changes  the  characteristics  of  the 
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X-ray  spectrum  until  it  no  longer  accurately  matches  the  assumed  22 
keV  exponential  spectrum.  The  ionization  and  heating  effects  become 
progressively  more  inaccurate  as  the  material  becomes  optically 
thick. 

This  presents  a  problem  when  dealing  with  the  wind  force 
shut-down  of  Model  3.  The  modified  spectrum  will  not  ionize  the 
oxygen  L-shells  under  the  same  conditions  assumed  for  our  model.  The 
50%  ionization  state  will  occur  closer  to  the  secondary  than  in  our 
model.  There  is  a  compensating  factor  which  must  also  be  taken  into 
account.  Hatchett,  Buff  and  McCray  (1976)  computed  the  ion 
populations  for  oxygen  in  an  optically  thick  gas.  As  noted 
previously,  the  oxygen  goes  from  neutral  to  fully  ionized  over  a  much 
smaller  range  of  £  than  in  the  optically  thin  case.  This  would  cause 
the  wind  force  to  turn  off  faster  and  more  closely  approximate  the 
sudden  turn  off  of  our  Model  3. 

Reconsidering  the  question  of  scaling,  recall  that  £  will  remain 
constant  as  the  mass  loss  rate  and  the  X-ray  luminosity  are  dropped. 
Decreasing  the  density  will  also  have  the  effect  of  decreasing  the 
optical  depths  in  the  model.  This  would  place  us  firmly  in  the  realm 
of  our  optically  thin  assumption.  We  conclude  that  Model  3  does  a 
fair  job  of  representing  the  physical  processes  which  occur  in  the 
real  system,  but  the  details  of  these  processes  are  limited  by  the 
assumptions  necessary  to  produce  the  model. 

E.  Comparisons  with  Previous  Works 

As  discussed  in  Chapter  I,  Livio,  Shara,  and  Shaviv  (1979) 
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(hereafter  refered  to  a  LSS)  presented  a  model  of  Cen  X-3  (4U1118-60) 
which  attempted  to  handle  the  hydrodynamics  and  X-ray  heating  in  a 
manner  similar  to  ours.  The  differences  between  the  two  models 
involves  the  method  in  which  the  X-ray  heating  is  accounted  for.  LSS 
assumed  a  constant  wind  velocity  with  v  =  v#  everywhere,  with  no 
provision  being  made  for  a  force  driving  the  wind.  Most  importantly, 
they  kept  the  X-ray  luminosity  fixed  instead  of  basing  it  on  the 
accretion  rate. 

LSS  cast  their  X-ray  heating  term  in  the  form  A(T-Te^),  where  A 
is  a  numerically  adjusted  parameter  serving  the  same  function  as  our 
t~*,  and  the  T-T  is  equivalent  to  our  E-E$s  (re.  Eq.  9.6).  The 
parameter  A  was  adjusted  numerically  until  thermal  balance  is  reached 
in  the  gas  and  is  then  allowed  to  remain  fixed.  It  seems  that  their 
method  of  fixing  the  parameter  A,  and  not  basing  the  X-ray  luminosity 
on  the  accretion  rate,  would  produce  an  artificially  smooth  problem 
with  no  means  of  supporting  the  thermal  instabilities  which  we 
observed.  For  their  heating  method  to  work,  they  had  to  assume  that 
the  flow  would  go  to  a  steady  state.  Their  model  could  be  compared 
with  our  Model  1  since  their  model  parameters  put  it  in  a  state  in 
which  little  effect  on  the  accretion  rate  due  to  X-ray  heating  would 
be  expected.  They  in  fact  predicted  no  instabilities  in  the  wind 
flows. 

The  velocity  plot  presented  by  LSS  shows  a  large  region  around 
the  primary.  It  seems  to  indicate  that  the  wind  has  almost  no 
velocity  in  the  region  perpendicular  to  the  line-of-centers  near  the 
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secondary.  The  velocity  of  the  wind  should  not  be  affected  far  from 
the  secondary  if  radiation  effects  on  the  wind  force  are  not 
considered.  In  fact,  the  velocity  at  the  side  boundaries  should  be 
the  unperturbed  velocity  of  the  wind  at  the  secondary's  orbit.  It 
appears  that  their  boundary  conditions  are  not  properly  supporting 
the  wind  velocities. 

Ignoring  the  wind  velocity  structure  also  adversely  impacts  on 
the  qualitive  features  of  their  results.  The  density  contours 
produced  by  LSS  do  not  show  the  pronounced  wake  which  we  found. 
Their  temperature  contours  should  be  nearly  spherical  in  front  of  the 
shock,  but  in  fact  are  very  irregular.  This  may  be  related  as  much 
to  their  method  of  X-ray  heating  as  to  their  wind  structure. 

Another  difference  between  the  two  models  is  the  treatment  used 
at  the  boundary  of  the  secondary.  We  accreted  material  through  the 
boundary  while  LSS  replaced  the  secondary  with  a  solid  boundary. 
This  caused  the  flow  to  build  up  to  hydrostatic  equilibrium  about  the 
secondary.  LSS  argued  that  this  approximates  the  case  where  the  mass 
flow  through  the  magnetosphere  of  the  secondary  is  slow  enough  so  as 
not  to  perturb  the  flow  in  the  immediate  vicinity.  It  is  not  clear 
how  this  slow  mass  flux  onto  the  secondary  produces  the  observed 
luminosity  and  still  does  not  perturb  the  incoming  flow.  To  support 
the  luminosity,  the  material  must  flow  to  the  surface  of  the  neutron 
star  at  essentially  the  same  rate  that  it  is  accreted. 

We  can  also  compare  our  results  to  those  found  by  Hoffman 
(1979).  Hoffman's  isothermal  model  developed  variations  in  the  mass 
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accretion  rate.  These  variations  were  a  result  of  modulating  the 
mass  flow  down  the  accretion  column  in  the  wake  of  the  secondary. 
Our  Models  1  and  2  did  not  develop  such  flow.  In  fact,  noticable 
accretion  columns  did  not  form  in  these  models.  Instead,  we  found 
the  bulk  of  the  accretion  flow  to  enter  from  the  sides  of  the  wake. 
Hoffmah's  models  did  not  include  the  stellar  wind  driving  force. 
This  additional  force  appears  to  be  the  reason  for  the  decrease  in 
the  accretion  flow  directly  behind  the  secondary  as  discussed  in 
Chapter  V.  His  models  also  did  not  show  the  formation  of  a  hot 
central  region  in  the  core  of  the  wake. 

Carlberg  (1978)  did  an  analytical  study  of  accretion  based  on 
the  line  accretion  model  which  included  arguments  concerning  the 
effects  of  X-ray  heating.  On  a  graph  of  the  wind  velocity  versus 
wind  density,  he  showed  regions  where  X-ray  heating  would  and  would 
not  affect  the  accretion  flow  (re.  his  Fig.  3.).  Since  he  assumed  a 
constant  wind  velocity,  we  can  assume  the  velocity  and  the  density  to 
be  the  wind  velocity  and  density  at  the  orbit  of  the  secondary.  On 
plotting  the  three  points  corresponding  to  our  models,  we  found  that 
all  three  models  lay  very  close  to  the  boundary  between  the  two 
regions.  Model  1  lay  just  marginally  below  the  boundary  in  the 
region  of  low  X-ray  impact.  Models  2  and  3  lay  just  above  the 
boundary,  predicting  significant  X-ray  impact  on  the  flows. 
Comparison  of  Model  3  is  difficult  in  light  of  the  variations  set  up 
in  the  wind.  However,  Models  1  and  2  appear  to  agree  with  his 
predictions.  The  drastic  change  in  accretion  rates  between  Models  1 


and  2  does  not  appear  to  be  justified  by  their  very  slight  difference 
in  position  relative  to  the  boundary.  This  may  again  be  the  result 
of  adding  the  wind  force  term  to  the  model. 

F.  Summary 

This  study  of  X-ray  heating  in  models  of  accretion  from  a 
stellar  wind  in  massive  X-ray  binary  stars  has  led  to  several 
conclusions.  From  Models  1  and  2  we  learned  that  the  analytic  line 
accretion  model  of  Hoyle  and  Lyttleton  (1939)  comes  within  a  factor 
of  10  of  describing  the  accretion  rates  in  these  systems.  The  degree 
to  which  the  line  accretion  model  agrees  with  the  X-ray  heated  models 
was  shown  to  be  dependent  on  the  X-ray  heating  rates  in  the  gas.  If 
the  X-ray  heating  is  sufficient  to  dominate  the  hydrodynamics  of  the 
flow,  then  the  accretion  rate  will  be  significantly  lower  than  that 
predicted  by  the  analytic  model. 

Our  ballistic  analysis  in  Chapter  VII  suggested  that  the  wind 
accretion  model  could  function  even  under  extreme  impairment  of  the 
wind  force  due  to  X-ray  ionization.  The  wind  is  severely  disrupted 
from  the  normal  spherical  outflow.  A  large  range  of  material 
velocities  was  found.  The  wind  flows  normally  outward  within  the 
X-ray  shadow  region.  Once  exposed  to  the  X-rays,  the  wind  will 
free-fall  under  the  influence  of  the  primary's  gravity.  Lower 
velocity  regions  will  begin  to  flow  back  in  towards  the  primary.  The 
analysis  suggests  that  sufficient  material  will  be  resent  in  the 
secondary's  vicinity  and  at  low  enough  velocity  to  allow  for 
efficient  accretion,  thus  easily  supporting  the  observed  luminosity. 


Material  flowing  inward,  however,  severely  distorts  the  spherically 
symmetric  outflow  assumed  by  the  radiation  driven  wind  model.  The 
impact  on  the  wind  model  is  not  known  at  this  time. 

In  Model  3  we  found  that  a  striking  feedback  mechanism  was 
created  in  a  flow  whose  wind  force  could  be  modified  by  the  X-ray 
radiation.  We  found  that  as  the  wind  was  turned  on  and  off, 
variations  in  density  were  created  which  caused  fluctuations  in  the 
X-ray  luminosity.  The  fluctuating  luminosity  led  to  variations  in 
the  wind  density.  The  characteristics  of  the  resultant  X-ray  flaring 
were  very  similiar  to  the  slow  X-ray  flares  observed  in  41)1700-37. 

These  results  pave  the  way  for  future  work  in  this  area.  A 
method  is  needed  to  simulate  three  dimensional  effects  efficiently 
and  accurately.  A  more  detailed  radiation  treatment  is  also 
necessary.  This  would  permit  a  more  accurate  determination  of  the 
X-ray  impact  on  the  wind,  and  would  make  it  possible  to  more 
accurately  predict  light  curves  for  the  binary  systems.  Such 
light-curves  would  provide  a  closer  observational  check  of  the 
models. 
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